MLAB Modeling of Interleukin13 Binding Kinetics
The MLAB advanced mathematical and statistical modeling system is a useful tool
for mathematical modeling as is demonstrated by the following example. MLAB is a
generalpurpose system which can be applied to model enzyme kinetics, multiple
site binding equilibrium, or any of a wide variety of other situations. Many
such examples are given at the website http://www.civilized.com/.
In this paper, we review some of MLAB's capabilities in a practical context of
independent interest.
Spaces between the cells that makeup the human body are filled with molecules that instigate a myriad of actions and reactions that form the biochemical activity of life. Some of these molecules, such as calcium and potassium ions, flow through channels within cell membranes to control muscle and heart cell contraction, or to act as the stuff of nerve impulses. Other molecules, such as estrogen and testosterone, have a more global role in many physiological subsystems of the body.
Mathematical modeling, either gross compartmental modeling, or specific physical chemistry modeling, such as chemical ligandbinding equilibrium or kinetic modeling, is one of the most powerful research tools available, either for quantifying the reactions involving specific ``signaling'' molecules or pharmaceutical agents that bind to receptor sites (which often occur as molecules embedded in cell membranes) or for helping to elucidate the plausible biochemical reactions that might be occurring. An example of this latter use of modeling to explore biochemical mechanism is discussed below.
Interleukin13(IL13) is one of the most important of the many signaling molecules that operate as a part of the immune system. IL13 binds to receptors on the surface of Blymphocytes, causing them to produce immunoglobulinG molecules. Moreover, IL13 binds to receptors on many other types of cells, including nonimmunesystem cells, to help produce the extremelycomplex socalled inflammation response. IL13 also binds to receptors on certain tumor cells, and, for a time, both controls their proliferation and instigates their death.
The antitumor effect of IL13 has been studied by many researchers [1],[2]. A recent study by Vladimir Kuznetsov, Nicholas Obiri, and Raj Puri, in which mathematical modeling played a central role, has shed light on IL13 binding with receptors on cells of a particular line of renal carcinoma [3]. This study involved a number of aspects of IL13 behavior interacting with a variety of receptors, but here we focus only on the binding of IL13 to receptors on the membranes of one particular type of tumor cell.
The data gathered specified the concentration of IL13 bound to receptors on the cellmembranes of a fixed known mass of tumor cells at various times during the binding reaction. After due care to convert cpm units to picomole(pM) units, subtract the effects of nonspecific binding, and deal with other complicating effects, we have the following data for four kinetic binding experiments, labeled experiment 1, experiment 2, experiment 3 and experiment 4. For experiment 1, 15 pM of IL13 ligand was added to .5×10^{6} tumor cells containing an unknown concentration of receptors. For experiment 2, 70 pM of IL13 ligand was added to .5×10^{6} tumor cells, for experiment 3, 200 pM of IL13 ligand was added to .5×10^{6} tumor cells, and for experiment 4, 500 pM of IL13 ligand was added to .5×10^{6} tumor cells. Thus, the same unknown concentration of receptors is assumed to be present in each experiment. In each experiment, the concentration in picomoles of bound IL13 ligand was measured at various times, as listed below.
time exp.1 exp.2 exp.3 exp.4 _________________________________ 1:      2: 0.03   14  3: 0.05  4.9   4: 0.06 1.5    5: 0.08    75 6: 0.17   36  7: 0.18  12.6   8: 0.25 2.66    9: 0.32    105 0: 0.33  16.1   11: 0.4   48  12: 0.41 3.12    13: 0.47    120 14: 0.56 3.6    15: 0.58   50  16: 0.67  17.5   17: 0.73 3.6    18: 0.75   54 135 19: 1 3.6 18.9  150 20: 1.25   58 155 21: 1.5 3.75    22: 1.67  20.3  145 23: 2   52  24: 3  18.2   25: 3.5 3.75    26: 5 3.75  52  27: 5.25  18.2   28: 6.2    150 29: 9   52 
Also, three socalled dissociation experiments, labeled experiment 5, experiment 6, and experiment 7, were done, where we started with a known amount of IL13 bound in various states, and then, in effect, removed it as it dissociated, and observed the concentration of IL13 measured in picomoles that remained bound at various times. This data was:
time exp.5 exp.6 exp.7 _______________________________________ 1: 0 3.96 52.8 157 2: 0.17  44.9328  3: 0.25 3.23928  138.0344 4: 0.33  46.3056  5: 0.42   139.416 6: 0.5 3.19176   7: 0.62   137.061 8: 0.67  44.4048  9: 0.75 3.1284   10: 1  44.1936 132.6179 11: 1.1 3.1284   12: 1.4   130.31 13: 1.7 3.12048   14: 2  43.2432 129.839 15: 3   129.996 16: 4  40.6032  17: 4.17 2.90664   18: 5 2.772  123.245 19: 6   120.262 20: 7   114.1861 21: 7.5  37.3296  22: 9   111.156 23: 11.8  34.32  24: 12.5 2.49084  
The MLAB mathematical and statistical modeling system was used to check a variety of possible models by fitting them to the data. It is important to note that, in each case, all seven sets of data were fit simultaneously using a combined set of model functions with shared parameters These models included one, two, and three independent receptor sites, as well as many other plausible scenarios.
Each model was fit to the data by estimating the unknown parameters in the model as those values that minimized the sum of the squared deviations of the model predictions from the observed data values. This process itself required careful attention to the mathematical formulation of the model, to the weightvalues employed, to the initial parameter guesses used, and to the constraints imposed on the parameter values to help guide the fitting process. The models often involved stiff differential equations which required the use of MLAB for their effective treatment.
Since neither the simple onesite model or twosite model adequately fit the data, a cooperative binding model was tried. In that model, it was postulated that the affinity constant characterizing the binding of the IL13 ligand to a receptor was a (linear) function of the concentration of the bound ligand, so that the binding affinity decreased (or increased) over time while binding was occurring. Fitting this model showed cooperativity appeared to be present and produced a better fit. In an attempt to explain this cooperativity effect, various other possibilities were explored. The result of weeks of work is presented below.
Kuznetsov endedup proposing that a third molecule, dubbed the coreceptor exists, and moreover, that the binding of IL13 with both the primary receptor and the coreceptor to form a ternary complex caused the release of a steady stream of instances of a fourth molecule called the inhibitor. These inhibitor molecules drift away from the cell membrane and reversibly bind to free IL13 molecules (this binding serves to inactivate the IL13 molecules that are thus captured so they cannot bind to the receptor to mediate any further associated biochemical events.) Since the inhibitor+IL13 ligand complex is assumed to be unable to bind to primary receptor molecules, the IL13 binding is effectively selfregulating; such binding provokes the appearance of inhibitor, which in turn, inhibits further binding of IL13 to the primary receptor! (This idea is consistent with the fact that larger doses of IL13 appear to have almost the same binding effect as smaller doses. It also implies that the inhibitor (if it exists) may need to be blocked in order to use IL13 as a therapeutic agent.)
It is important to note that the postulated molecules can be replaced by other mechanisms. For example, the coreceptor could, in fact, be merely a conformational change (possibly mediated in some manner) which allows a reshaped bound complex molecule itself to link to either another IL13 molecule or, perhaps, a free primary receptor molecule so as to inactivate its target. The inhibitor may be produced via an enzyme reaction where the IL13+receptor+coreceptor ternary complex acts as an enzyme, or, perhaps, the ternary complex acts to open a ``channel'' in the cell membrane that allows the release of inhibitor molecules.
Although the modeling reported here has been fruitful in allowing some mechanisms to be rejected as implausible, and in showing that certain other mechanisms may be involved, further biochemical analysis is required to make a more convincing case for any such specific suggestions. For example, it is important to note that the model functions might be ``over fit'', due to the levels of uncertainty in the observed data, even though the visual results are impressive. The role of modeling here has thus been to narrow the plausible kinds of mechanisms that might be sought. We cannot expect more definitive results without more data from experiments tailored to expose specific mechanisms, and the modeling presented here helps us design such experiments. Note, however, that in the case where we are modeling such a biochemical schema in order to quantify some particular behavior, like a clearance rate or a metabolic conversion rate, the exact mechanism may be less important than the accurate quantification provided by a suitable model. (Models for glucoseinsulin interaction fall in this category; they are physiologically crude, but practically useful.)
Let L denote the IL13 ligand, let R denote the primary receptor, and let C denote the coreceptor. Also, let B denote the L+R complex, let G denote the B+C complex, let H denote the inhibitor released due to the action of the G material, and let X denote the (inhibited) L+H complex.
Thus, we have Kuznetsov's model:
L + R « B
B + C « G
G ® G + H
L + H « X
In both the kinetic binding experiments and the dissociation experiments, the
concentration of the total bound IL13 ligand was observed, i.e.
the combined concentration of B and G. In fact, there still appears to be a
small amount of cooperativity present in the L +
R « B reaction. This is modeled by
assuming that the association constant for this reaction is of the form k_{1}
L^{b} with b < 1, which decreases as the concentration of free
IL13 decreases during the time of reaction. Note, however the initial
association constant increases as the initial amount of IL13 ligand
supplied increases. Kuznetsov explains this cooperativity effect as an effect
due to the rough nature of the cell membrane which imposes a slowing of binding
due to the difficulty of diffusing the ligand into the ``crevices'' of the cell
membrane. (This diffusion effect is only apparent when several sets of kinetic
experimental data are fit simultaneously; a single set of data can be fit by a
variety of simpler models, but the schema proposed here is the only one found
that accomodates all the data.) Also, note that this model is surely not
complete; there is no mechanism postulated to scavenge the inhibitor for
example, so that, as stated, our system would eventually be flooded with
inhibitor, and virtually all the IL13 ligand would be bound with inhibitor.
Let L_{j}(t) denote the concentration of free IL13 ligand at
time t in the kinetic binding experiment j and let L0_{j} = L_{j}(0),
the initial concentration of free IL13 ligand. For our four binding
experiments, we have L0_{1}=15, L0_{2}=70, L0_{3}=200,
and L0_{4}=500.
Let R_{j}(t) denote the concentration of free primary receptor at
time t in the kinetic binding experiment j and let R0 denote the initial
concentration of free primary receptor, which we take to be identical in all our
experiments.
Let C_{j}(t) denote the concentration of free coreceptor at time
t in the kinetic binding experiment j and let C0 denote the initial
concentration of free coreceptor, which we take to be identical in all our
experiments.
Let H_{j}(t) denote the concentration of free inhibitor at time t
in the kinetic binding experiment j. We assume that H_{j}(0)=0.
Let B_{j}(t) denote the concentration of IL13+receptor bound complex at
time t in the kinetic binding experiment j. We assume that B_{j}(0)=0.
Let G_{j}(t) denote the concentration of IL13+receptor+coreceptor
bound complex at time t in the kinetic binding experiment j. We assume that G_{j}(0)=0.
Let X_{j}(t) denote the concentration of IL13+inhibitor bound complex
at time t in the kinetic binding experiment j. We assume that X_{j}(0)=0.
Finally, let A_{j}(t) denote the ratio of the concentration of membrane
bound IL13 ligand in either the form B or the form G at time t to the total
concentration of IL13 ligand in the kinetic binding experiment j, so that A_{j}(t)=(B_{j}(t)+G_{j}(t))/L0_{j}
with A_{j}(0)=0.
We use the fraction bound because it is easier to compare the binding curves
from different experiments using fraction units. We have four such kinetic
binding experiments, so j=1,2,3,4.
Then, following the above stochiometric reactions, the algebraic and differential equations that define these functions are:

The quantities k_{1} L_{j}^{b} and k_{1} are the association and dissociation ``constants'' for the reaction L+ R « B. The quantities k_{2} and k_{2} are the association and dissociation constants for the reaction B + C « G. The quantity k_{3} is the association constant for the reaction G ® G + H. The quantities k_{4} and k_{4} are the association and dissociation constants for the reaction L + H « X.
The values of L0_{1}, L0_{2}, L0_{3}, and L0_{4} are the known values 15, 70, 200, and 500, respectively. The parameters to be estimated are: R0, C0, k_{1}, b, k_{1}, k_{2}, k_{2}, k_{3}, k_{4}, k_{4}.
The three dissociation experiments are described by:
B ® L + R
B + C « G
Again, for these dissociation experiments, the concentration of the total bound IL13 ligand was observed, i.e. the combined concentration of B and G, was measured at various times while continually ``washing'' away the nonmembrane bound species.
Recall that we label the three dissociation experiments as experiments number 5,6, and 7. Then for j=5,6,7, we will have sets of functions that compose our models for the three sets of dissociation data.
In fact, dissociation experiment 5 is based on kinetic binding experiment 1; in particular, dissociation experiment 5 started when a quasisteadystate was reached after about 5 hours of binding occurred in a copy of kinetic binding experiment 1. That means we started with a concentration of L0_{1} micromoles of IL13 ligand, R0 micromoles of receptor, C0 micromoles of coreceptor, and 0 micromoles of inhibitor, and after t_{d}:=5 hours, we began washing off the nonmembrane bound material and measuring the bound IL13 ligand. Thus the initial concentrations of B and G in experiment 5 are whatever was ``inherited'' from experiment 1! We know the sum of these concentrations is 3.96 pM, since a measurement was made at the time of starting the dissociation experiment, but we do not know the individual concentrations of B and G separately.
However, since MLAB allows auxilliary functions defined in terms of other differentialequationdefined functions to appear in a set of functions being fit to data, We can estimate these component concentrations separately by fitting the function ZB_{5}(t)=B_{1}(t)B0_{5} to the point (t_{d}, 0) and the function ZG_{5}(t)=G_{1}(t)G0_{5} to the point (t_{d}, 0) by varying the parameters B0_{5} and G0_{5}! The parameters B0_{5} and G0_{5} will be thereby computed as estimates of the values of B_{1}(t_{d}) and G_{1}(t_{d}), and these estimates may be used as the initial values for the timecourse functions that estimate the concentrations of B and G in the model equations for dissociation experiment 5. Note the constraint B0_{5}+G0_{5}=3.96 should be approximately correct.
Dissociation experiment 6 is based on kinetic binding experiment 3; in particular, as with experiment 5, dissociation experiment 6 started after about 5 hours of binding occurred in a copy of kinetic binding experiment 3. That means we started with a concentration of L0_{3} micromoles of IL13 ligand, R0 micromoles of receptor, C0 micromoles of coreceptor, and 0 micromoles of inhibitor, and after reaching quasisteadystate near time t_{d}, we began washing off the nonmembrane bound material and measuring the bound IL13 ligand. Thus the initial concentrations of B and G in experiment 6 are whatever was ``inherited'' from experiment 3. The sum of these initial concentrations was 52.8 pM.
Again, we can estimate these component concentrations separately by fitting the function ZB_{6}(t)=B_{3}(t)B0_{6} to the point (t_{d}, 0) and the function ZG_{6}(t)=G_{3}(t)G0_{6} to the point (t_{d}, 0) by varying the parameters B0_{6} and G0_{6}. The parameters B0_{6} and G0_{6} will be thereby computed as estimates of the values of B_{3}(t_{d}) and G_{3}(t_{d}), and these estimates may be used as the initial values for the timecourse functions that estimate the concentrations of B and G in the model equations for dissociation experiment 6. Note the constraint B0_{6}+G0_{6}=52.8 should be approximately correct.
Finally, dissociation experiment 7 is based on kinetic binding experiment 4; again, dissociation experiment 7 started after about 5 hours of binding occurred in a copy of kinetic binding experiment 7; thus we started with a concentration of L0_{4} micromoles of IL13 ligand, R0 micromoles of receptor, C0 micromoles of coreceptor, and 0 micromoles of inhibitor, and after about t_{d} hours, we began washing off the nonmembrane bound material and measuring the bound IL13 ligand. Thus the initial concentrations of B and G in experiment 7 are whatever was ``inherited'' from experiment 4. The sum of these initial concentrations was 157 pM.
Just as before, we can estimate these component concentrations separately by fitting the function ZB_{7}(t)=B_{4}(t)B0_{7} to the point (t_{d}, 0) and the function ZG_{7}(t)=G_{4}(t)G0_{7} to the point (t_{d}, 0) by varying the parameters B0_{7} and G0_{7}. The parameters B0_{7} and G0_{7} will be thereby computed as estimates of the values of B_{4}(t_{d}) and G_{4}(t_{d}), and these estimates may be used as the initial values for the timecourse functions that estimate the concentrations of B and G in the model equations for dissociation experiment 7. Note the constraint B0_{7}+G0_{7}=157 should be approximately correct.
Now for the three dissociation experiments numbered 5,6, and 7, we may take
j=5,6,7, and state the following functions that compose our model for the
dissociation data.
Let B_{j}(t) denote the concentration of IL13+receptor bound complex at
time t in the kinetic dissociation experiment j. We specify the initial value B_{j}(0)
to be B0_{j}, which will be estimated as discussed above.
Let G_{j}(t) denote the concentration of IL13+receptor+coreceptor
bound complex at time t in the kinetic dissociation experiment j. We specify the
initial value G_{j}(0) to be G0_{j}, which will be estimated as
discussed above.
Let L_{j}(t) denote the concentration of free IL13 ligand at
time t in the kinetic dissociation experiment j. We assume that L_{j}(0)=0.
Let C_{j}(t) denote the concentration of free coreceptor at time
t in the kinetic dissociation experiment j Thus C_{j}(t) = C0G_{j}(t),
where C0 is the inital concentration of coreceptor (before any IL13 ligand is
added.)
Finally, let A_{j}(t) denote the ratio of the concentration of membrane
bound IL13 ligand in either the form B or the form G at time t to the total
concentration of IL13 ligand in the kinetic dissociation experiment j, so that
A_{j}(t)=(B_{j}(t)+G_{j}(t))/(B0_{j}+G0_{j}).
Then the algebraic and differential equations that define these functions are:

The quantity k_{1} is the dissociation constant for the reaction B ® L + R. The quantities k_{2} and k_{2} are the association and dissociation constants for the reaction B + C « G. Note the constants C0, k_{1}, k_{2}, and k_{2} also appear in the models for the kinetic binding experiments. The additional parameters B0_{5}, B0_{6}, B0_{7}, G0_{5}, G0_{6}, and G0_{7} also appear in these functions as initial conditions. They are defined by the auxillary functions:

The parameters to be estimated are: C0, k_{1}, k_{2}, k_{2}, B0_{5}, B0_{6}, B0_{7}, G0_{5}, G0_{6}, and G0_{7}.
Overall then, we have seven sets of data to be fit by functions defined by seven corresponding sets of differentialalgebraic equations, with 16 unknown shared parameters C0, R0, b, k_{1}, k_{1}, k_{2}, k_{2}, k_{3}, k_{4}, k_{4}, B0_{5}, B0_{6}, B0_{7}, G0_{5}, G0_{6}, and G0_{7} to be estimated. The MLAB commands (which were entered in a dofile) and the results are given below. It is important to note that this is a synthesized presentation; the actual modeling required many dofiles, and many strategic steps to obtain the final results. Also, please note the small mathematical devices introduced below. In particular, we define the function p to avoid computing the nonintegral power of a (small) negative number; in theory, L_{j} should never be less than 0, but numerical error can be introduced during solving that violates this constraint. Also, note the use of appropriate weights and constraints in the MLAB dofile below.
/* dofile: il13.do = fit Kuznetsov's model to experimental data */ fct p(L,b)=if L<0 then 0 else L^b /* Model functions for Experiment 1 */ fct B1't(t) = k1*p(L1(t),1+b)*(R0B1G1) km1*B1 k2*B1*(C0G1) +km2*G1 fct G1't(t) = k2*B1*(C0G1) km2*G1 fct H1't(t) = k3*G1 k4*L1(t)*H1 +km4*X1 fct X1't(t) = k4*L1(t)*H1 km4*X1 fct L1(t) = L01B1G1X1 fct A1(t) = (B1+G1)/L01 init B1(0)=0 init G1(0)=0 init H1(0)=0 init X1(0)=0 /* Model functions for Experiment 2 */ fct B2't(t) = k1*p(L2(t),1+b)*(R0B2G2) km1*B2 k2*B2*(C0G2) +km2*G2 fct G2't(t) = k2*B2*(C0G2) km2*G2 fct H2't(t) = k3*G2 k4*L2(t)*H2 +km4*X2 fct X2't(t) = k4*L2(t)*H2 km4*X2 fct L2(t) = L02B2G2X2 fct A2(t) = (B2+G2)/L02 init B2(0)=0 init G2(0)=0 init H2(0)=0 init X2(0)=0 /* Model functions for Experiment 3 */ fct B3't(t) = k1*p(L3(t),1+b)*(R0B3G3) km1*B3 k2*B3*(C0G3) +km2*G3 fct G3't(t) = k2*B3*(C0G3) km2*G3 fct H3't(t) = k3*G3 k4*L3(t)*H3 +km4*X3 fct X3't(t) = k4*L3(t)*H3 km4*X3 fct L3(t) = L03B3G3X3 fct A3(t) = (B3+G3)/L03 init B3(0)=0 init G3(0)=0 init H3(0)=0 init X3(0)=0 /* Model functions for Experiment 4 */ fct B4't(t) = k1*p(L4(t),1+b)*(R0B4G4) km1*B4 k2*B4*(C0G4) +km2*G4 fct G4't(t) = k2*B4*(C0G4) km2*G4 fct H4't(t) = k3*G4 k4*L4(t)*H4 +km4*X4 fct X4't(t) = k4*L4(t)*H4 km4*X4 fct L4(t) = L04B4G4X4 fct A4(t) = (B4+G4)/L04 init B4(0)=0 init G4(0)=0 init H4(0)=0 init X4(0)=0 /* Model functions for Experiment 5 */ fct B5't(t) = km1*B5 k2*B5*(C0G5) +km2*G5 fct G5't(t) = k2*B5*(C0G5) km2*G5 fct A5(t) = (B5+G5)/(B05+G05) fct ZB5(t) = B1(t)B05 fct ZG5(t) = G1(t)G05 init B5(0)=B05 init G5(0)=G05 /* Model functions for Experiment 6 */ fct B6't(t) = km1*B6 k2*B6*(C0G6) +km2*G6 fct G6't(t) = k2*B6*(C0G6) km2*G6 fct A6(t) = (B6+G6)/(B06+G06) fct ZB6(t) = B3(t)B06 fct ZG6(t) = G3(t)G06 init B6(0)=B06 init G6(0)=G06 /* Model functions for Experiment 7 */ fct B7't(t) = km1*B7 k2*B7*(C0G7) +km2*G7 fct G7't(t) = k2*B7*(C0G7) km2*G7 fct A7(t) = (B7+G7)/(B07+G07) fct ZB7(t) = B4(t)B07 fct ZG7(t) = G4(t)G07 init B7(0)=B07 init G7(0)=G07 /* Set constants */ L01=15; L02=70; L03=200; L04=500 /* Set parameter guesses near Kuznetsov's best values */ R0= 2990 C0= 2990 B= 0.109050928 K1= 0.0003001336855 KM1= 4.8 K2= 0.0001216742301 KM2= 0.026 K3= 0.7673768901 K4= 5 KM4= 0.07631566379 B05= 0.8712 G05= 3.0888 B06= 10.15 G06= 47.85 B07= 21 G07= 129 /* Readin the data and construct the data matrices D1D7, and Z */ bd=read(bind,100,5); dd=read(diss,100,4); D1=0 compress(bd col (1,2),2); D1 col 2=(D1 col 2)/L01 D2=0 compress(bd col (1,3),2); D2 col 2=(D2 col 2)/L02 D3=0 compress(bd col (1,4),2); D3 col 2=(D3 col 2)/L03 D4=0 compress(bd col (1,5),2); D4 col 2=(D4 col 2)/L04 D5=compress(dd col (1,2),2); D5 col 2=(D5 col 2)/D5[1,2] D6=compress(dd col (1,3),2); D6 col 2=(D6 col 2)/D6[1,2] D7=compress(dd col (1,4),2); D7 col 2=(D7 col 2)/D7[1,2] delete bd,dd Z=5&'0 /* Define the constraint set Q */ constraints Q={ \ B05+G05>3.5, B05+G05<4.5, \ B06+G06>44, B06+G06<65, \ B07+G07>140, B07+G07<170, \ b>0, b<.13, \ R0>1500, R0<3250, \ C0>1500, C0<3250, \ k1>0, k1<.01, \ km1>3, km1<5, \ k2>0, \ km2>.02,km2<.03, \ k3>.1, k3<2, \ k4>.5, k4<5, km4>0, km4<5 } /*Setup weight vectors */ wt1=d1 col 2; wt1[1]=wt1[6] wt2=d2 col 2; wt2[1]=wt2[6] wt3=d3 col 2; wt3[1]=wt3[6] wt4=d4 col 2; wt4[1]=wt4[6] wz5=list(1/3.96) wz6=list(1.33/58) wz7=list(2/150) /* Do the curvefitting  fit our model to the weighted data. */ errfac=.00001 method=gear maxiter=40 lsqrpt =11 fit(R0,C0,b,k1,km1,k2,km2,k3,k4,km4,B05,G05,B06,G06,B07,G07), \ A1 to D1 with wt wt1, \ A2 to D2 with wt wt2, \ A3 to D3 with wt wt3, \ A4 to D4 with wt wt4, \ A5 to D5, \ A6 to D6, \ A7 to D7, \ ZB5 to Z with wt wz5, \ ZG5 to Z with wt wz5, \ ZB6 to Z with wt wz6, \ ZG6 to Z with wt wz6, \ ZB7 to Z with wt wz7, \ ZG7 to Z with wt wz7, constraints Q
The output of this dofile produced by the fit statement follows.
Begin iteration 1 bestsosq=6.006547e001 MarquardtLevenberg line search probe: eps = 1.000000e010 sum of squares = 0.04582580707113 Begin iteration 2 bestsosq=4.582581e002 MarquardtLevenberg line search probe: eps = 1.000000e010 sum of squares = 0.01204906773736 Begin iteration 3 bestsosq=1.204907e002 MarquardtLevenberg line search probe: eps = 1.000000e010 sum of squares = 0.01204619061608 MarquardtLevenberg line search probe: eps = 1.624368e003 sum of squares = 0.01204625000369 MarquardtLevenberg line search probe: eps = 1.624368e001 sum of squares = 0.01204635763616 MarquardtLevenberg line search probe: eps = 1.624368e+001 sum of squares = 0.01204821700603 MarquardtLevenberg line search probe: eps = 1.624368e+002 sum of squares = 0.01204894914491 MarquardtLevenberg line search probe: eps = 1.624368e+003 sum of squares = 0.0120490553849 final parameter values value error dependency parameter 3250 12182.85387 0.9999997438 R0 3250 13462.31276 0.9999999084 C0 0.1049589281 0.05167059669 0.999952039 B 0.000296334322 0.0010658965 0.999999698 K1 5 0.6762822057 0.9997546225 KM1 0.000106107999 0.00044342822 0.999999997 K2 0.02302618885 0.001314873569 0.9332693779 KM2 0.814427359 0.08304489826 0.9996683989 K3 5 229.3106734 0.9998679795 K4 0.3069022341 30.49902014 0.9998620474 KM4 0.8770250533 0.06970645033 0.90829221013 B05 3.144308218 0.2154997399 0.983383752 G05 9.141654224 0.7304766809 0.9838603169 B06 46.6525215 2.285711761 0.9982697999 G06 21.15065076 1.60517971 0.9940830812 B07 118.8493492 6.223070522 0.9995980458 G07 3 iterations; CONVERGED best weighted sum of squares = 1.204619e002 weighted root mean square error = 1.441156e002 weighted deviation fraction = 1.429964e002 lagrange multiplier[5] = 4.110273827e006 lagrange multiplier[10] = 5.569797706e008 lagrange multiplier[12] = 7.365789198e008 lagrange multiplier[16] = 0.001165085076 lagrange multiplier[23] = 2.823709132e006
The graphic results of the above curvefitting are shown below. MLAB was used to produce the pictures presented below, but we have elided the actual commands involved. The starting parameters above are Kuznetsov's final result; here we have adjusted all sixteen parameters simultaneously with different weights for the data than those used by Kuznetsov to obtain a different fit. The fact that many different sets of parameter values can be found which yield good fits under the same or slightly different conditions is to be expected. This modeling problem is illconditioned, and, in particular, there are high correlations among certain of the parameters that makes the determination of the ``true values'' of those parameters impossible. Other parameters, however, do not change very much, i.e. are stable across various fitting experiments, and we have somewhat more faith in these parameters. The important point here is that Kuznetsov's model does admit a good fit with physicallyplausible parameter values, whereas other competing models failed to do so. We can find other, mathematically''better'' fits if we relax our constraints, but these fits have implausible parameter values.
0
The picture below shows the nature of the kinetic curves of Kuznetsov's model in the particular case of 200 pM of IL13 initially introduced. The solidline curve is B, the dashedline curve is G, and the dottedline curve is X. (Do you think the fact that the three kinetic curves presented appear to intersect at the same point has any significance?)
References
[1] Debinski, W.; Obiri, N.I.; Pastan,I.; Puri, R.K., ``A Novel Chimeric Protein Composed of Interleukin13 and Pseudomonas exotoxin is highly cytotoxic to Human Carcinoma Cells Expressing Receptors for Interleukin4'', J. Biol. Chem. 270, pp. 1677516780, 1995
[2] Obiri, N.I.; Husain ,S.R.,; Debinski, W.; Puri, R.K., ``Interleukin13 Inhibits Growth of Renal Cell Carcinoma Cells Independently of the p140 Interleukin4 Receptor Chain'', Clinical Cancer Res. 2, pp. 17431749, 1996
[3] Kuznetsov, Vladimir A.; Puri, Raj K., ``Kinetic Analysis of HighAffinity Forms of Interleukin13 Receptors: Suppression of IL13 Binding by IL2 Receptor g Chain'', Biophysical Journal, Vol. 77, pp. 154172, 1999