Fitting multiple parameters to individual data sets while simultaneously fitting their combination to another data set

I have 4 parameters (A-D) that combine into another parameter Z as follows:

    Z = A - (B + C - B.*C.*D/1336)                                          %(1)

I also have ample measured data for A(t), B(t), C(t) and Z(t), but unfortunately most available datasets only provide one or two parameters simutaneously, almost never all 5.

Individually, parameters A-C can well be fitted well their respective data sets using these functions:

    Model_A = a1*exp(a2*t) + a3*exp(a4*t) + a0;                             %(2a)
    Model_B = b1*exp(b2*t) + b3*exp(b4*t) + b0;                             %(2b)
    Model_C = c1*exp(c2*t) + c3*exp(c4*t) + c0;                             %(2c)

And for D there is a theoretical function of the form:

    Model_D = d1*exp(d2*t) + d0;                                            %(2d)

For parameter Z there does not seem to be such a neat function. Only a 6th order polynomial seems to do a reasonable job, but I don't like it as it lacks physical meaning.

Inputing fits (2a-d) into equation (1) seems to get the general shape of Z right, but, as can be expected, this is far from optimized (see figure).

So that is why I was wondering whether it would be possible to somehow link models (2a-d) in such a way that they would each be optimized for their individual data sets, while simultaneaously optimizing for the data of parameter Z.

Is there a way to accomplish this in Matlab?