numerics - How to fit 3 data sets to a model of 4 differential equations?

I'm a biologist and a newbie in Mathematica. I want to fit three data sets to a model consisting of four differential equations and 10 parameters. I want to find the parameters best fitting to my model. I have searched the forum and found several related examples. However, I could not find anything that matched my question.

Here are the details:

I have three time-series datasets: (xdata, ydata, zdata)

time = Quantity[{0, 3, 7, 11, 18, 25, 38, 59}, "seconds"];
tend = QuantityMagnitude[Last[time]];

xdata:

xdata = Quantity[{0, 0.223522, 0.0393934, 0.200991, 0.786874, 1, 
    0.265464, 0.106174}, "milligram"];

xfitdata = QuantityMagnitude[Transpose[{time, xdata}]];

ydata:

ydata = Quantity[{0, 0.143397, 0.615163, 0.628621, 0.53515, 0.519805, 
    0.757092, 1}, "milligram"];
yfitdata = QuantityMagnitude[Transpose[{time, ydata}]];

wdata:

wdata = Quantity[{0.0064948, 0.221541, 1, 0.434413, 0.732392, 
    0.458638, 0.1484432, 0.0294298}, "milligram"];

wfitdata = QuantityMagnitude[Transpose[{time, wdata}]];

I used ParametricNDSolve to solve the 4-DE model:

pfun = {x, y, z, w} /.  
  ParametricNDSolve[{x'[t] == 
     k1 - k10 x[t] w[t - 25] - k2 x[t] - k3 w[t] w[t], 
    y'[t] == -k8 y[t] + k10 x[t] w[t - 25] + k3 w[t] x[t], 
    z'[t] == k4 y[t] - k5 z[t], 
    w'[t] == (k6 x[t])/(y[t]^n + 1) - k7 w[t], x[t /; t <= 0] == 0.01,
     y[t /; t <= 0] == 0.01, z[t /; t <= 0] == 0.01, 

    w[t /; t <= 0] == 0.01}, {x, y, z, w}, {t, 0, tend}, {k1, k2, k3, 
    k4, k5, k6, k7, k8, n, k10}]

Then I used FindFit. But I don't know how to specify that xdata is supposed to be fitted to x[t], zdata to z[t] and wdata to w[t] via least-squares fit. For y[t], there are no time-series data, but the parameter (k8) for y[t] is supposed to be determined as well.

I have tried the following, which is apparently wrong:

fit = FindFit[xfitdata, 
  pfun[{k1, k2, k3, k4, k5, k6, k7, k8, n, k10}][
   t], {{k1, 0.0859}, {k2, 0.0125}, {k3, 0.8541}, {k4, 0.0185}, {k5, 
    0.1004}, {k6, 0.5002}, {k7, 0.0511}, {k8, 0.0334}, {n, 9}, {k10, 
    0.8017}}, t]

This is the error message:

FindFit::nrlnum: The function value {0. +<<1>>[0.],-0.223522+<<1>>,-0.0393934+<<1>>,-0.200991+<<1>>,-0.786874+<<1>>[{0.0859,0.0125,0.8541,0.0185,0.1004,0.5002,0.0511,0.0334,9.,0.8017}][18.],-1.+<<1>>[25.],-0.265464+<<1>>,-0.106174+<<1>>[59.]} is not a list of real numbers with dimensions {8} at {k1,k2,k3,k4,k5,k6,k7,k8,n,k10} = {0.0859,0.0125,0.8541,0.0185,0.1004,0.5002,0.0511,0.0334,9.,0.8017}. >>

I'm lost and I would really appreciate your help!

Answer

Since the question isn't clear about which datasets are which and arguably has too many parameters, I'll use the example from here instead:

$$ \begin{array}{l} A+B\underset{k_2}{\overset{k_1}{\leftrightharpoons }}X \\ X+B\overset{k_3}{\longrightarrow }\text{products} \\ \end{array} \Bigg\} \Longrightarrow A+2B\longrightarrow \text{products} $$

We solve the system and generate some fake data:

sol = ParametricNDSolveValue[{

    a'[t] == -k1 a[t] b[t] + k2 x[t], a[0] == 1,
    b'[t] == -k1 a[t] b[t] + k2 x[t] - k3 b[t] x[t], b[0] == 1,
    x'[t] == k1 a[t] b[t] - k2 x[t] - k3 b[t] x[t], x[0] == 0
    }, {a, b, x}, {t, 0, 10}, {k1, k2, k3}
   ];

abscissae = Range[0., 10., 0.1];
ordinates = With[{k1 = 0.85, k2 = 0.15, k3 = 0.50},
   Through[sol[k1, k2, k3][abscissae], List]
  ];


data = ordinates + RandomVariate[NormalDistribution[0, 0.1^2], Dimensions[ordinates]];
ListLinePlot[data, DataRange -> {0, 10}, PlotRange -> All, AxesOrigin -> {0, 0}]

The data look like this, where blue is A, purple is B, and gold is X:

The key to the exercise, of course, is the simultaneous fitting of all three datasets in order for the rate constants to be determined self-consistently. To achieve this we have to prepend to each point a number, i, that labels the dataset:

transformedData = {
    ConstantArray[Range@Length[ordinates], Length[abscissae]] // Transpose,
    ConstantArray[abscissae, Length[ordinates]],

    data
   } ~Flatten~ {{2, 3}, {1}};

We also need a model that returns the values for either A, B, or X depending on the value of i:

model[k1_, k2_, k3_][i_, t_] := 
  Through[sol[k1, k2, k3][t], List][[i]] /;
    And @@ NumericQ /@ {k1, k2, k3, i, t};

The fitting is now straightforward. Although it will help if reasonable initial values are given, this is not strictly necessary here:

fit = NonlinearModelFit[

   transformedData,
   model[k1, k2, k3][i, t],
   {k1, k2, k3}, {i, t}
  ];

The result is correct. Worth noting, however, is that the off-diagonal elements of the correlation matrix are quite large:

fit["CorrelationMatrix"]
(* -> {{ 1.,        0.764364, -0.101037},
       { 0.764364,  1.,       -0.376295},

       {-0.101037, -0.376295,  1.      }} *)

Just to be sure of having directly addressed the question, I will note that the process does not change if we have less than the complete dataset available (although the parameters might be determined with reduced accuracy in this case). Typically it will be most difficult experimentally to measure the intermediate, so let's get rid of the dataset for X (i == 3) and try again:

reducedData = DeleteCases[transformedData, {3, __}];
fit2 = NonlinearModelFit[
   reducedData,
   model[k1, k2, k3][i, t],
   {k1, k2, k3}, {i, t}
  ];

The main consequence is that the error on $k_3$ is significantly larger:

This can be seen to be the result of greater correlation between $k_1$ and $k_3$ when fewer data are available for fitting:

fit2["CorrelationMatrix"]
(* -> {{ 1.,        0.7390200,  -0.1949590},
       { 0.7390200, 1.,          0.0435416},
       {-0.1949590, 0.0435416,   1.       }} *)

On the other hand, the correlation between $k_2$ and $k_3$ is greatly reduced, so that all of the rate constants are still sufficiently well determined and the overall result does not change substantially.