How can I perform sensitivity analysis on ODEs with SALib?

Question

def iron(XYZ,t,a12,a21,a23,a32,b13,b31,I):
X1,X2,X3=XYZ
dX1=-a12*(X1)+a21*(X2)-b13*(X1)+b31*(X3)
dX2=-a23*(X2)-a21*(X2)+a12*(X1)+a32*(X3)
dX3=-a32*(X3)-b31*(X3)+a23*(X2)+b13*(X1)-I
return dX1,dX2,dX3;
a12=0.0005 
a21=0.00001
a23=0.0003 
a32=0.0002 
b13=0.0001 
b31=0.000001 
I=0.001 

XYZ0=[1000.,30.,10.]
X10=1000.
X20=50.
X30=30.

t=linspace(0,100,1000) #(start,stop,num samples to generate)
XYZ=odeint(iron,XYZ0,t,args=(a12,a21,a23,a32,b13,b31,I))

Is it possible to perform sensitivity analysis on this system of ODEs using SALib? I want to study the influence of model inputs (parameters a and b, initial condition). Also, can I obtain the asymptotic solution values?

Answer 1

The following code is a sample implementation of the problem for Sobol analysis. The bounds for each parameter would have to be adjusted based on the problem; I assumed a range for the example. Provided there is an asymptotic solution that can be found generally as a function of the parameter (a12, a21,...,b13,b31,I), you could follow a similar procedure. Determining the asymptotic solution is probably better posted as a separate question.

The X1, X2, and X3 values at the end of the time must be analyzed separately by the Sobol method. The sensitivity index could be computed for each X1, X2, and X3 for all time steps, but it would require saving all the output from each iteration of the loop. It would also require running the Sobol analysis many more times.

Some of the sample output for the following code example is:

  ====X2 Sobol output====

  Parameter S1 S1_conf ST ST_conf
    a12 0.409635 0.083264 0.411180 0.049683
    a21 0.000002 0.000095 0.000001 0.000000
    a23 -0.000626 0.002955 0.000471 0.000057
    a32 0.000068 0.000504 0.000017 0.000002
    b13 0.000045 0.000232 0.000004 0.000001
    b31 0.000000 0.000000 0.000000 0.000000
    x1_0 0.430008 0.078269 0.434074 0.053487
    x2_0 0.169098 0.051591 0.162944 0.018678
    x3_0 -0.000038 0.000335 0.000007 0.000001

Example Code Implementation

# importing packages
from scipy import integrate as sp
import numpy as np
import SALib
from SALib.sample import saltelli
from SALib.analyze import sobol

# definition of the system of ODEs
def iron(XYZ,t,a12,a21,a23,a32,b13,b31,I):
  X1,X2,X3=XYZ
  dX1=-a12*(X1)+a21*(X2)-b13*(X1)+b31*(X3)
  dX2=-a23*(X2)-a21*(X2)+a12*(X1)+a32*(X3)
  dX3=-a32*(X3)-b31*(X3)+a23*(X2)+b13*(X1)-I
  return dX1,dX2,dX3;

# default parameter values
a12=0.0005 
a21=0.00001
a23=0.0003 
a32=0.0002 
b13=0.0001 
b31=0.000001 
I=0.001 

# initial condition
XYZ0=[1000.,30.,10.]
X10=1000.
X20=50.
X30=30.

# tmie steps
t=np.linspace(0,100,1000) #(start,stop,num samples to generate)

# example single calculation
XYZ=sp.odeint(iron,XYZ0,t,args=(a12,a21,a23,a32,b13,b31,I))

### Sobol analysis ###
# defining problem
# can add the 'I' parameter
# assumed that range for each parameter is 80-120% of value assumed above
# can be changed
problem = {
  'num_vars': 9, #a's, b's and initial condition
  'names': ['a12', 'a21','a23','a32','b13','b31','x1_0','x2_0','x3_0'],
  'bounds':  np.column_stack((np.array([a12, a21,a23,a32,b13,b31,XYZ0[0],XYZ0[1],XYZ0[2]])*0.8,np.array([a12, a21,a23,a32,b13,b31,XYZ0[0],XYZ0[1],XYZ0[2]])*1.2))
}

# Generate samples
vals = saltelli.sample(problem, 500)

# initializing matrix to store output
Y = np.zeros([len(vals),1])

# Run model (example)
# numerically soves the ODE
# output is X1, X2, and X3 at the end time step
# could save output for all time steps if desired, but requires more memory
Y = np.zeros([len(vals),3])
for i in range(len(vals)):
  Y[i][:] = sp.odeint(iron,[vals[i][6],vals[i][7],vals[i][8]],t,\
    args=(vals[i][0],vals[i][1],vals[i][2],vals[i][3],vals[i][4],vals[i][5],I))[len(XYZ)-1]


# completing soboal analysis for each X1, X2, and X3
print('\n\n====X1 Sobol output====\n\n')
Si_X1 = sobol.analyze(problem, Y[:,0], print_to_console=True)
print('\n\n====X2 Sobol output====\n\n')
Si_X2 = sobol.analyze(problem, Y[:,1], print_to_console=True)
print('\n\n====X3 Sobol output====\n\n')
Si_X3 = sobol.analyze(problem, Y[:,2], print_to_console=True)