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

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.

Note Added in proof

Oleksander's answer provides a better fit to the data than my solution below and circumvents the reiteration-and-solving-individually problem I describe below.

This is not an answer to your question specifically, but rather is one method to use ParametricNDSolve to fit experimental data.

Defining the problem

$\require{mhchem}$ Differential equations come in to play in chemical kinetics, most notably in determining the rates of reactions. Consider the following reaction $\ce{A -> B + C}$. To my knowledge, the set of equations that describe the differential rate laws for the three components of this system cannot be solved analytically. Let's say that we have collected a set of concentration vs. time data for each of the three species A, B, and C and would like to fit the data to this chemical equation. First, I will make some data, noting that if the rate constant for the reverse reaction is k then the rate of the forward reaction can be given by K x k where K is the equilibrium constant for the chemical equation.

tdata = NDSolve[{a'[t] == -b'[t] == -c'[t] == -k K a[t] + k b[t] c[t],
     a[0] == a0, b[0] == b0, c[0] == c0} /. {K -> 3, k -> 1, 
    a0 -> 0.5, b0 -> 0.4, c0 -> 0.1}, {a[t], b[t], c[t]}, {t, 0, 1}]
edata = Flatten[
   Table[{t, RandomReal[{0.98, 1.02}] a[t], 
      RandomReal[{0.98, 1.02}] b[t], 
      RandomReal[{0.98, 1.02}] c[t]} /. tdata, {t, 0, 1, 0.05}], 1];

Here, I've chosen the answers to be k = 1, K = 3, and the initial concentrations of A, B, and C to be 0.5, 0.4, and 0.1, respectively.

Solving the problem

We start by using ParametricNDSolve to generate the interpolation functions, and generate a plot just to make sure we are on the right track.

sol = ParametricNDSolve[{a'[t] == -b'[t] == -c'[t] == -k K a[t] + 
     k b[t] c[t], a[0] == a0, b[0] == b0, c[0] == c0}, {a, b, c}, {t, 
   0, 1}, {k, K, a0, b0, c0}];
a1 = a[1, 3, 0.5, 0.4, 0.1] /. sol;
b1 = b[1, 3, 0.5, 0.4, 0.1] /. sol;
c1 = c[1, 3, 0.5, 0.4, 0.1] /. sol;
Plot[Evaluate@{a1[t], b1[t], c1[t]}, {t, 0, 1}]

Using NonlinearModelFit with ParametricNDsolve solutions

The solutions from ParametricNDSolve can be used directly in the NonlinearModelFit function, providing results that seem to provide a decent fit. Note that fitting the concentration data for B and C threw errors (more on that below).

nlma = NonlinearModelFit[edata[[All, {1, 2}]], 
  a[k, K, a0, b0, c0][t] /. sol, {k, K, a0, b0, c0}, t]
nlmb = NonlinearModelFit[edata[[All, {1, 3}]], 
  b[k, K, a0, b0, c0][t] /. sol, {k, K, a0, b0, c0}, t]
nlmc = NonlinearModelFit[edata[[All, {1, 4}]], 
  c[k, K, a0, b0, c0][t] /. sol, {k, K, a0, b0, c0}, t]

The results look pretty good, however the parameters are actually meaningless:

Adding constraints

This problem can be improved by adding constraints to the NonlinearModelFit

nlma2 = NonlinearModelFit[
  edata[[All, {1, 2}]], {a[k, K, a0, b0, c0][t] /. sol2, k > 0, K > 0,
    0 < a0 < 1, 0 < b0 < 1, 0 < c0 < 1}, {k, K, a0, b0, c0}, t]
nlmb2 = NonlinearModelFit[
  edata[[All, {1, 3}]], {b[k, K, a0, b0, c0][t] /. sol2, k > 0, K > 0,
    0 < a0 < 1, 0 < b0 < 1, 0 < c0 < 1}, {k, K, a0, b0, c0}, t]
nlmc2 = NonlinearModelFit[
  edata[[All, {1, 4}]], {c[k, K, a0, b0, c0][t] /. sol2, k > 0, K > 0,
    0 < a0 < 1, 0 < b0 < 1, 0 < c0 < 1}, {k, K, a0, b0, c0}, t]

The fit is better, but still not great, and it is important to note that the data sets for each individual concentration give varying optimal parameters. I'm a fan of the brute force and ignorance approach, and one way to address the problem with the parameters is to take the Mean and StandardDeviation of the "optimal" parameters and use these as new constraints.

(* New stuff *)
mean = Mean[{k, K, a0, b0, c0} /. #["BestFitParameters"] & /@ {nlma2, 
    nlmb2, nlmc2}]
sd = StandardDeviation[{k, K, a0, b0, c0} /. #[
      "BestFitParameters"] & /@ {nlma2, nlmb2, nlmc2}]
const = MapThread[{#1 - #2 < #3 < #1 + #2} &, {mean, 
   sd, {k, K, a0, b0, c0}}]
(* old stuff *)
nlma3 = NonlinearModelFit[
  edata[[All, {1, 2}]], {a[k, K, a0, b0, c0][t] /. sol2, const}, {k, 
   K, a0, b0, c0}, t]
nlmb3 = NonlinearModelFit[
  edata[[All, {1, 3}]], {b[k, K, a0, b0, c0][t] /. sol2, const}, {k, 
   K, a0, b0, c0}, t]
nlmc3 = NonlinearModelFit[
  edata[[All, {1, 4}]], {c[k, K, a0, b0, c0][t] /. sol2, const}, {k, 
   K, a0, b0, c0}, t]
TableForm[{k, K, a0, b0, c0} /. #["BestFitParameters"] & /@ {nlma3, 
   nlmb3, nlmc3}, 
 TableHeadings -> {{"a", "b", "c"}, {"k(1)", "K(3)", "a0(0.5)", 
    "b0(0.4)", "c0(0.1)"}}]

We are now starting to get agreement between the three sets of data. After repeating this loop approximately a dozen times, the results start to settle:

Conclusion

We can see that some of the best fit parameters are reasonably close to the actual values (concentrations of A and B) what might be considered the important values (k and K) have appreciable, but possibly acceptable errors (on the order of 20 to 30%). I presume, although have not tested, that part of the problem lies in redefining the constraints for concentration of C, which from the first iteration did not contain the true value. The conclusion to draw from this exercise is that a complex model (and this is not a terribly complex model) will give you a fit of the data, but you need to know something about the reality of the parameters before you can actually trust the fit. Maybe Mathematica version 10 will provide that insight, but as of now, we must use our own brains to determine the value of a fit.