# Deducing parameters of 2-by-2 linear system of first-order ODEs from data

I'm currently researching glucose-insulin regulation system and I came across this homogeneous system of linear differential equations: $$\begin{array}{lr} x'&=&-\alpha x-\beta y\\ y'&=&\gamma x-\delta y \end{array}$$ Here $$x$$ and $$y$$ are the difference in the concentrations of glucose and insulin from their equilibrium levels respectively. $$\alpha$$ is the rate constant related to the efficiency that the liver absorbs glucose, $$\beta$$ to the rate at which glucose is absorbed by muscle, $$\gamma$$ to the rate that insulin is produced by the pancreas and $$\delta$$ the rate at which insulin is degraded by the liver.

I want to find the parameters $$\alpha, \beta, \gamma,$$ $$\delta$$, but I have absolutely no idea how to do so. I tried to integrate them and plug in the initial conditions and other values, but I wasn't able to find them. Can you please help me find these constants? Thank you.

Note: the data is summarized in the table below

P.S. In the case that this problem cannot be solved by hand, would it be possible to use computer to solve it?

• If we let $v=\begin{pmatrix} x \\ y\end{pmatrix}$ and $A=\begin{pmatrix} -\alpha & -\beta \\ \gamma & -\delta\end{pmatrix}$, then this system may be written as $v'=Av$ with formal solution $v(t)=e^{At}v(0)$. The trouble is that the matrix exponential $e^{At}$ is in general a quite complicated function of time, so this doesn't exactly solve the problem in a practical sense. Commented Jan 20, 2023 at 17:44
• What do you mean by more definition? Currently, the only information that I have is every parameters have to be greater than 0. Commented Jan 20, 2023 at 19:09
• For reference, the closed form solution (i.e. the afore-mentioned matrix exponential) is given here (using $a,b,c,d$ instead of $\alpha,\beta,\gamma,\delta$). (That said I don't follow the objection being made.) Commented Jan 20, 2023 at 20:56
• Typically done numerically. In practice you have six parameters ($\alpha$, $\beta$, $\gamma$, $\delta$, $x_0$, $y_0$) which are unknowns. For any choice of these, compute the solution at the measured times (for example simple Runge-Kutta integration) and compare the computed results with the measured, say via mean squared error which then provides a basis for optimizing the parameters. This approach works for more complex, nonlinear systems also. Commented Jan 20, 2023 at 23:44
• Cross-post from stackoverflow.com/questions/75188506/…. Please only ask once or at least link the questions. I gave a list of SO topics of fitting ODE parameters using scipy.odr and scipy.optimize. There are frameworks like gekko or pyomo that are more systematic and general than ODR. Commented Jan 21, 2023 at 10:27

Follows a MATHEMATICA script which shows the fitting process to obtain the parameters $$(\alpha, \beta,\gamma,\delta)$$ and initial conditions $$(c_1,c_2)$$. After determining the solution closed form, the code is as plain as possible, being illustrative of each step.

Clear[X, Y, ERROR, alpha, beta, gamma, delta, c1, c2, d, d2]
tk = {0.0, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 4, 6};
y1 = {70, 150, 165, 145, 90, 75, 65, 75, 80, 75};
y2 = {8, 60, 44, 32, 15, 9, 8, 6, 2, 1};
sol = DSolve[{x'[t] == -alpha x[t] - beta y[t],
y'[t] == gamma x[t] - delta y[t]}, {x, y}, t][[1]] /. {C[1] -> c1, C[2] -> c2};
{xt, yt} = {x[t], y[t]} /. sol // FullSimplify

X[t_?NumberQ, alpha_?NumberQ, beta_?NumberQ, gamma_?NumberQ, delta_?NumberQ, c1_?NumberQ, c2_?NumberQ] := Module[{d2, d},
d2 = (alpha - delta)^2 - 4 beta gamma;
If[d2 < 0, d2 = -d2];
d = Sqrt[d2];
Return[(E^(-(1/2) (alpha + d + delta) t) (2 beta c2 + c1 (alpha + d - delta) + (-2 beta c2 + c1 (-alpha + d + delta)) E^(d t)))/(2 d)]
]

Y[t_?NumberQ, alpha_?NumberQ, beta_?NumberQ, gamma_?NumberQ, delta_?NumberQ, c1_?NumberQ, c2_?NumberQ] := Module[{d2, d},
d2 = (alpha - delta)^2 - 4 beta gamma;
If[d2 < 0, d2 = -d2]; d = Sqrt[d2];
Return[(E^(-(1/2) (alpha + d + delta) t) (c2 (-alpha + d + delta) -  2 c1 gamma + E^(d t) (c2 (alpha + d - delta) + 2 c1 gamma)))/(2 d)]
]

ERROR[alpha_?NumberQ, beta_?NumberQ, gamma_?NumberQ, delta_?NumberQ, c1_?NumberQ, c2_?NumberQ] :=
Module[{},
Return[Sum[Abs[X[tk[[k]], alpha, beta, gamma, delta, c1, c2] - y1[[k]]]/(2 + tk[[k]]) + Abs[Y[tk[[k]], alpha, beta, gamma, delta, c1, c2] - y2[[k]]]/(2 + tk[[k]]), {k, 1, 10}]]]

vars = {alpha, beta, gamma, delta, c1, c2};
sol2 = NMinimize[{ERROR[alpha, beta, gamma, delta, c1, c2], alpha > 0,beta > 0, gamma > 0, delta > 0, (alpha - delta)^2 - 4 beta gamma > 10^-9}, vars, Method -> "DifferentialEvolution"]
{xt, yt} = {x[t], y[t]} /. sol /. sol2[[2]] // N;

gr1 = Plot[{xt, yt}, {t, 0, 6}];
gr2 = ListPlot[Transpose[{tk, y1}]];
gr3 = ListPlot[Transpose[{tk, y2}]];
Show[gr1, gr2, gr3]


NOTE

The involved functions

$$\cases{ X(t,\alpha,\beta,\gamma,\delta,c_1,c_2) = \frac{e^{-\frac{1}{2} t (\alpha +d+\delta )} \left(e^{d t} (c_1 (-\alpha +d+\delta )-2 \beta c2)+c_1 (\alpha +d-\delta )+2 \beta c_2\right)}{2 d}\\ Y(t,\alpha,\beta,\gamma,\delta,c_1,c_2) =\frac{e^{-\frac{1}{2} t (\alpha +d+\delta )} \left(e^{d t} (c_2 (\alpha +d-\delta )+2 c_1 \gamma )+c_2 (-\alpha +d+\delta )-2c_1 \gamma \right)}{2 d}}$$

with $$d = \sqrt{(\alpha -\delta )^2-4 \beta \gamma }$$ and

$$ERROR(\alpha,\beta,\gamma,\delta,c_1,c_2) = \sum_k\left( |X(t_k,\alpha,\beta,\gamma,\delta,c_1,c_2)-x_k| + |Y(t_k,\alpha,\beta,\gamma,\delta,c_1,c_2)-y_k|\right)/(2+t_k)$$

• Your excellent plot reveals an interesting pattern: that one of the terms being modeled appears to approach a nonzero equilibrium value. In practice this means that two additional parameters need to be included in the model ( since x and y are defined as deviations from these two unknown equilibria). Commented Jan 21, 2023 at 0:45
• @Cesareo: very nice script (+1).
– Moo
Commented Jan 23, 2023 at 15:38

Off hand I would suggest that you (i) use your table of numerical data to try to estimate the time rates of change that appear on the left side of your model such as $$\frac{dx}{dt}\approx \frac{ \Delta x}{\Delta t}$$; then (ii) attempt to fit a pair of linear equations in $$x$$ and $$y$$ to that estimated data. The standard method for performing (ii) is called linear regression (based on least-squared error minimization). You could collaborate with someone versed in statistics to execute the details.

• This will work if the time step can be made sufficiently small. Unfortunately, the time steps in the data set given seem to be rather large and thus not likely to be well-approximated by linear fits. Commented Jan 20, 2023 at 20:57