Reverse problem of finding time-varying parameters of an ODE with the help of solution data.

Any way in which I can find the time-varying parameters of an ODE ; The solution to which is know in the form of time-series data.

Background:

Example: SIR model ode: as following.

% The SIR model differential equations.

def deriv(y, t, beta, gamma):

S, I, R = y

dSdt = -beta * S * I

dIdt = beta * S * I - gamma * I

dRdt = gamma * I

return dSdt, dIdt, dRdt

Then, I use odeint(deriv, y0, t, args=( beta, gamma)) (In python) to solve for S, I and R using minimisation of these parameters.

Question

But I want to ask, is there any way in which a reverse problem can be constructed such that; I have the data for S, I and R as time-sereis. and then we can calculate 'beta' and 'gamma' paramter from there?

Note beta and gamma should be time dependent.They can be assumed as sum of natural cubic splines; Where N_k(t), k=1 to K, are K natural cubic spline basis functions evaluated at K-2 equally spaced knots in addition to the boundary knots .

Can this problem be solved on matlab?

Thanks.

[Torsten]

Aren't S, I and R time-dependent enough to keep the parameters beta and gamma fix ?

[Anshuman]

Yes they are time dependent as I have written in deriv function. But I want to make changes by letting beta and gamma to be time dependent as well! As in the phenomenon they are time dependent.

Also, I wish to assume that log(β(t)) = ∑K natural cubic splines. (I saw this in a research paper) ; Is there any way I can formulate this problem mathematically?

[Torsten]

Is there any way I can formulate this problem mathematically?

Maybe, but for turning beta and gamma in time-dependent functions, you will need many many experimental data for S, I and R, I guess.

[Anshuman]

I have time series data set for S, I and R values. A data for [S,I,R] for 100 time-steps; I now wish to estimate the time varying parameters beta and gamma.

Maybe a mre simpler assumption of beta = a1*(t) + b1 and gamma = a2*(t) + b2. and now estimate a1,b1,a2,b2 ? Or the above mentioned K-spline method for beta and delta.

I am not able to devise a code to formulate all this.

[Anshuman]

I have written a code for estimating 'beta' and ;delta' values for the whole time step; that is it will return a constant beta and delta value for the whole timespan in python.

Code:

(a) Define the ODE model of SIR with cosntant parameter 'beta' and 'gamma'

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

I0, R0 = 0.10798003802281400, 0.0020386203150461700
S0 = 1 - I0 - R0
t = np.linspace(0, 146, 146)

# The SIR model differential equations.

def deriv(y, t, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I 
    dIdt = beta * S * I  - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt

I_d = x[:,1] # from dataset, we will finally try to fit the SIR plot with this data 

(b) Define fn(x) for minimising the loss with respect to the data of 'I'

def fn(x):
# parameters unwrapped
    beta, gamma = x;   
    I0, R0 = 0.10798003802281400, 0.0020386203150461700
    S0 = 1 - I0 - R0
    y0 = S0, I0, R0
    ret = odeint(deriv, y0, t_d, args=(beta, gamma))
    S, I, R = ret.T
    error = sum((I-I_d)*(I-I_d))

    return error

# initialise with current best guess
init_x = [0.1, 0.05]

res = minimize(fn, init_x, method='Nelder-Mead', tol=1e-8)

fn(res.x) # calculate final error
res.x # calculate final parameters

(c) Finally we obtain a constant 'beta' and delta' values for the full time-span.

[Torsten]

Fine. And why do you want to change a winning python ?

[Anshuman]

Actually I couldn't find a way to implement python packages for this time dependent problem; so I thought of switching to matlab.

Answer 1

[Torsten]

Use "lsqcurvefit" to fit the parameters together with an integrator (e.g. ODE45).

For an example, see Star Strider's code under

https://de.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

Here, four parameters from a Monod Kinetic are fitted against measurement data.

[Anshuman]

Are these paramwts tie dependent?

[Torsten]

For your first shot, you wanted to estimate the constants a1, b1, a2 and b2, and these are not time-dependent. So the above code can be applied.

For fully time-dependent estimates, look this video to understand how to proceed:

https://www.youtube.com/watch?v=60xX4R708Iw

[Anshuman]

Sure, thanks a lot!