MLAB Modeling of Interleukin-13 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 general-purpose 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 web-site 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 make-up 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 sub-systems of the body.

Mathematical modeling, either gross compartmental modeling, or specific physical chemistry modeling, such as chemical ligand-binding 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.

Interleukin-13(IL-13) is one of the most important of the many signaling molecules that operate as a part of the immune system. IL-13 binds to receptors on the surface of B-lymphocytes, causing them to produce immunoglobulin-G molecules. Moreover, IL-13 binds to receptors on many other types of cells, including non-immune-system cells, to help produce the extremely-complex so-called inflammation response. IL-13 also binds to receptors on certain tumor cells, and, for a time, both controls their proliferation and instigates their death.

The anti-tumor effect of IL-13 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 IL-13 binding with receptors on cells of a particular line of renal carcinoma [3]. This study involved a number of aspects of IL-13 behavior interacting with a variety of receptors, but here we focus only on the binding of IL-13 to receptors on the membranes of one particular type of tumor cell.

The data gathered specified the concentration of IL-13 bound to receptors on the cell-membranes 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 non-specific 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 IL-13 ligand was added to .5106 tumor cells containing an unknown concentration of receptors. For experiment 2, 70 pM of IL-13 ligand was added to .5106 tumor cells, for experiment 3, 200 pM of IL-13 ligand was added to .5106 tumor cells, and for experiment 4, 500 pM of IL-13 ligand was added to .5106 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 IL-13 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 so-called dissociation experiments, labeled experiment 5, experiment 6, and experiment 7, were done, where we started with a known amount of IL-13 bound in various states, and then, in effect, removed it as it dissociated, and observed the concentration of IL-13 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 weight-values 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 one-site model or two-site 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 IL-13 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 ended-up proposing that a third molecule, dubbed the co-receptor exists, and moreover, that the binding of IL-13 with both the primary receptor and the co-receptor 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 IL-13 molecules (this binding serves to inactivate the IL-13 molecules that are thus captured so they cannot bind to the receptor to mediate any further associated biochemical events.) Since the inhibitor+IL-13 ligand complex is assumed to be unable to bind to primary receptor molecules, the IL-13 binding is effectively self-regulating; such binding provokes the appearance of inhibitor, which in turn, inhibits further binding of IL-13 to the primary receptor! (This idea is consistent with the fact that larger doses of IL-13 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 IL-13 as a therapeutic agent.)

It is important to note that the postulated molecules can be replaced by other mechanisms. For example, the co-receptor 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 IL-13 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 IL-13+receptor+co-receptor 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 glucose-insulin interaction fall in this category; they are physiologically crude, but practically useful.)

Let L denote the IL-13 ligand, let R denote the primary receptor, and let C denote the co-receptor. 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 + H

L + H    X


In both the kinetic binding experiments and the dissociation experiments, the concentration of the total bound IL-13 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 k1 Lb with b < 1, which decreases as the concentration of free IL-13 decreases during the time of reaction. Note, however the initial association constant increases as the initial amount of IL-13 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 IL-13 ligand would be bound with inhibitor.

Let Lj(t) denote the concentration of free IL-13 ligand at time t in the kinetic binding experiment j and let L0j = Lj(0), the initial concentration of free IL-13 ligand. For our four binding experiments, we have L01=15, L02=70, L03=200, and L04=500.
Let Rj(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 Cj(t) denote the concentration of free co-receptor at time t in the kinetic binding experiment j and let C0 denote the initial concentration of free co-receptor, which we take to be identical in all our experiments.
Let Hj(t) denote the concentration of free inhibitor at time t in the kinetic binding experiment j. We assume that Hj(0)=0.
Let Bj(t) denote the concentration of IL-13+receptor bound complex at time t in the kinetic binding experiment j. We assume that Bj(0)=0.
Let Gj(t) denote the concentration of IL-13+receptor+co-receptor bound complex at time t in the kinetic binding experiment j. We assume that Gj(0)=0.
Let Xj(t) denote the concentration of IL-13+inhibitor bound complex at time t in the kinetic binding experiment j. We assume that Xj(0)=0.
Finally, let Aj(t) denote the ratio of the concentration of membrane bound IL-13 ligand in either the form B or the form G at time t to the total concentration of IL-13 ligand in the kinetic binding experiment j, so that Aj(t)=(Bj(t)+Gj(t))/L0j with Aj(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:


dBj
dt
(t)
=
(k1 Ljb)Lj Rj -k-1Bj - dGj
dt
dGj
dt
(t)
=
k2 Bj Cj -k-2Gj
dHj
dt
(t)
=
k3 Gj - dXj
dt
dXj
dt
(t)
=
k4 Lj Hj - k-4Xj
Lj(t)
=
L0j-Bj(t)-Gj(t)-Xj(t)
Rj(t)
=
R0-Bj(t)-Gj(t)
Cj(t)
=
C0-Gj(t)
Aj(t)
=
(Bj(t)+Gj(t))/L0j

The quantities k1 Ljb and k-1 are the association and dissociation ``constants'' for the reaction L+ R    B. The quantities k2 and k-2 are the association and dissociation constants for the reaction B + C    G. The quantity k3 is the association constant for the reaction  G + H. The quantities k4 and k-4 are the association and dissociation constants for the reaction L + H    X.

The values of L01, L02, L03, and L04 are the known values 15, 70, 200, and 500, respectively. The parameters to be estimated are: R0, C0, k1, b, k-1, k2, k-2, k3, k4, k-4.

The three dissociation experiments are described by:

  L + R

B + C    G

Again, for these dissociation experiments, the concentration of the total bound IL-13 ligand was observed, i.e. the combined concentration of B and G, was measured at various times while continually ``washing'' away the non-membrane 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 quasi-steady-state 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 L01 micromoles of IL-13 ligand, R0 micromoles of receptor, C0 micromoles of co-receptor, and 0 micromoles of inhibitor, and after td:=5 hours, we began washing off the non-membrane bound material and measuring the bound IL-13 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 differential-equation-defined functions to appear in a set of functions being fit to data, We can estimate these component concentrations separately by fitting the function ZB5(t)=B1(t)-B05 to the point (td, 0) and the function ZG5(t)=G1(t)-G05 to the point (td, 0) by varying the parameters B05 and G05! The parameters B05 and G05 will be thereby computed as estimates of the values of B1(td) and G1(td), and these estimates may be used as the initial values for the time-course functions that estimate the concentrations of B and G in the model equations for dissociation experiment 5. Note the constraint B05+G05=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 L03 micromoles of IL-13 ligand, R0 micromoles of receptor, C0 micromoles of co-receptor, and 0 micromoles of inhibitor, and after reaching quasi-steady-state near time td, we began washing off the non-membrane bound material and measuring the bound IL-13 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 ZB6(t)=B3(t)-B06 to the point (td, 0) and the function ZG6(t)=G3(t)-G06 to the point (td, 0) by varying the parameters B06 and G06. The parameters B06 and G06 will be thereby computed as estimates of the values of B3(td) and G3(td), and these estimates may be used as the initial values for the time-course functions that estimate the concentrations of B and G in the model equations for dissociation experiment 6. Note the constraint B06+G06=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 L04 micromoles of IL-13 ligand, R0 micromoles of receptor, C0 micromoles of co-receptor, and 0 micromoles of inhibitor, and after about td hours, we began washing off the non-membrane bound material and measuring the bound IL-13 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 ZB7(t)=B4(t)-B07 to the point (td, 0) and the function ZG7(t)=G4(t)-G07 to the point (td, 0) by varying the parameters B07 and G07. The parameters B07 and G07 will be thereby computed as estimates of the values of B4(td) and G4(td), and these estimates may be used as the initial values for the time-course functions that estimate the concentrations of B and G in the model equations for dissociation experiment 7. Note the constraint B07+G07=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 Bj(t) denote the concentration of IL-13+receptor bound complex at time t in the kinetic dissociation experiment j. We specify the initial value Bj(0) to be B0j, which will be estimated as discussed above.
Let Gj(t) denote the concentration of IL-13+receptor+co-receptor bound complex at time t in the kinetic dissociation experiment j. We specify the initial value Gj(0) to be G0j, which will be estimated as discussed above.
Let Lj(t) denote the concentration of free IL-13 ligand at time t in the kinetic dissociation experiment j. We assume that Lj(0)=0.
Let Cj(t) denote the concentration of free co-receptor at time t in the kinetic dissociation experiment j Thus Cj(t) = C0-Gj(t), where C0 is the inital concentration of co-receptor (before any IL-13 ligand is added.)
Finally, let Aj(t) denote the ratio of the concentration of membrane bound IL-13 ligand in either the form B or the form G at time t to the total concentration of IL-13 ligand in the kinetic dissociation experiment j, so that Aj(t)=(Bj(t)+Gj(t))/(B0j+G0j).

Then the algebraic and differential equations that define these functions are:


dBj
dt
(t)
=
-k-1Bj - dGj
dt
dGj
dt
(t)
=
k2 Bj Cj -k-2Gj
Cj(t)
=
C0-Gj(t)
Aj(t)
=
(Bj(t)+Gj(t))/(B0j+G0j)

The quantity k-1 is the dissociation constant for the reaction B L + R. The quantities k2 and k-2 are the association and dissociation constants for the reaction B + C    G. Note the constants C0, k-1, k2, and k-2 also appear in the models for the kinetic binding experiments. The additional parameters B05, B06, B07, G05, G06, and G07 also appear in these functions as initial conditions. They are defined by the auxillary functions:


ZB5(t)
=
B1(t)-B05
ZG5(t)
=
G1(t)-G05
ZB6(t)
=
B3(t)-B06
ZG6(t)
=
G3(t)-G06
ZB7(t)
=
B4(t)-B07
ZG7(t)
=
G4(t)-G07

The parameters to be estimated are: C0, k-1, k2, k-2, B05, B06, B07, G05, G06, and G07.

Overall then, we have seven sets of data to be fit by functions defined by seven corresponding sets of differential-algebraic equations, with 16 unknown shared parameters C0, R0, b, k1, k-1, k2, k-2, k3, k4, k-4, B05, B06, B07, G05, G06, and G07 to be estimated. The MLAB commands (which were entered in a do-file) and the results are given below. It is important to note that this is a synthesized presentation; the actual modeling required many do-files, 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 non-integral power of a (small) negative number; in theory, Lj 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 do-file below.

/* do-file: 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)*(R0-B1-G1) -km1*B1 -k2*B1*(C0-G1) +km2*G1  
fct G1't(t) = k2*B1*(C0-G1) -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) = L01-B1-G1-X1 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)*(R0-B2-G2) -km1*B2 -k2*B2*(C0-G2) +km2*G2  
fct G2't(t) = k2*B2*(C0-G2) -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) = L02-B2-G2-X2  
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)*(R0-B3-G3) -km1*B3 -k2*B3*(C0-G3) +km2*G3  
fct G3't(t) = k2*B3*(C0-G3) -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) = L03-B3-G3-X3  
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)*(R0-B4-G4) -km1*B4 -k2*B4*(C0-G4) +km2*G4  
fct G4't(t) = k2*B4*(C0-G4) -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) = L04-B4-G4-X4  
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*(C0-G5) +km2*G5  
fct G5't(t) = k2*B5*(C0-G5) -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*(C0-G6) +km2*G6  
fct G6't(t) = k2*B6*(C0-G6) -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*(C0-G7) +km2*G7  
fct G7't(t) = k2*B7*(C0-G7) -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 


/* Read-in the data and construct the data matrices D1-D7, 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 } 


/*Set-up 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 curve-fitting - 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 do-file produced by the fit statement follows.


Begin iteration 1 bestsosq=6.006547e-001 
  Marquardt-Levenberg line search probe: eps = 1.000000e-010 
    sum of squares = 0.04582580707113 
Begin iteration 2 bestsosq=4.582581e-002 
  Marquardt-Levenberg line search probe: eps = 1.000000e-010 
     sum of squares = 0.01204906773736 
Begin iteration 3 bestsosq=1.204907e-002 
  Marquardt-Levenberg line search probe: eps = 1.000000e-010 
    sum of squares = 0.01204619061608 
  Marquardt-Levenberg line search probe: eps = 1.624368e-003 
     sum of squares = 0.01204625000369 
  Marquardt-Levenberg line search probe: eps = 1.624368e-001 
    sum of squares = 0.01204635763616 
Marquardt-Levenberg line search probe: eps = 1.624368e+001 
    sum of squares = 0.01204821700603 
  Marquardt-Levenberg line search probe: eps = 1.624368e+002 
    sum of squares = 0.01204894914491 
  Marquardt-Levenberg 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.204619e-002 
weighted root mean square error = 1.441156e-002 
weighted deviation fraction = 1.429964e-002 
lagrange multiplier[5] = -4.110273827e-006 
lagrange multiplier[10] = 5.569797706e-008 
lagrange multiplier[12] = 7.365789198e-008 
lagrange multiplier[16] = 0.001165085076 
lagrange multiplier[23] = 2.823709132e-006
 

The graphic results of the above curve-fitting 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 ill-conditioned, 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 physically-plausible 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 IL-13 initially introduced. The solid-line curve is B, the dashed-line curve is G, and the dotted-line 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 Interleukin-13 and Pseudomonas exotoxin is highly cytotoxic to Human Carcinoma Cells Expressing Receptors for Interleukin-4'', J. Biol. Chem. 270, pp. 16775-16780, 1995

[2] Obiri, N.I.; Husain ,S.R.,; Debinski, W.; Puri, R.K., ``Interleukin-13 Inhibits Growth of Renal Cell Carcinoma Cells Independently of the p140 Interleukin-4 Receptor Chain'', Clinical Cancer Res. 2, pp. 1743-1749, 1996

[3] Kuznetsov, Vladimir A.; Puri, Raj K., ``Kinetic Analysis of High-Affinity Forms of Interleukin-13 Receptors: Suppression of IL-13 Binding by IL-2 Receptor g Chain'', Biophysical Journal, Vol. 77, pp. 154-172, 1999