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

Question

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?

Answer 1

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.

Answer 2

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) $$