A MATHEMATICAL MODEL FOR FIBROBLAST GROWTH FACTOR COMPETITION BASED ON ENZYME KINETICS

In this paper, we develop a mathematical model for the competition of two species of fibroblast growth factor, FGF-1 and FGF-2, for the same cell surface receptor. We provide pathways for this interaction using experimental data obtained by Neufeld and Gospodarowicz reported in 1986 [9]. These pathways demonstrate how the interaction of two fibroblast growth factors affects cell proliferation. Upon development of these pathways, we use simulations in MATLAB and optimization to extrapolate the values of a variety of biochemical parameters imbedded within the model. Furthermore, it should be possible to use the model as the basis for a testable hypothesis. We explore this predictive ability with further simulations in MATLAB.

1. Introduction.Fibroblast growth factors (FGFs), among the earliest growth factors to be identified and purified, constitute a large family of at least 25 unique but related secreted proteins that stimulate cell proliferation and that are expressed in many tissues.The levels of FGFs found in these tissues are regulated by many biological factors, which reflects the involvement of the FGFs in more than one physiological event.One important physiological function of FGFs is in wound healing.The role of FGFs in wound healing has been demonstrated in many ways, among them FGFs are found in wound fluids, the absence of FGF2 delays wound healing [6], and the expression of many FGF genes increases after wounding [8].In combination with other growth factors, FGFs also play many roles in early embryonic development including to define the dorso-ventral pattern of the neural tube [14], to promote limb development [11], and to define the structure of the early embryo [3].The FGFs act through specific receptors (FGFR) that initiate signals inside the cell to alter cellular functions such as gene expression.The importance of these receptors to normal development is demonstrated by the many human skeletal diseases caused by mutations in FGFR genes [15].
2. Biological Activity of FGF-1 and FGF-2.In [9], Neufeld and Gospodarowicz examined the physical and chemical characteristics of FGF-1 and FGF-2. 1 Noting the apparent similarities between FGF-1 and FGF-2, Neufeld and Gospodarowicz proceeded to investigate the differential affinities of these two FGFs to the same cell surface receptor.Several experiments were carried out to characterize biological activity of FGF-1 and FGF-2 in a variety of different situations.Experiment 1. First, the effects of increasing concentrations of either FGF-1 and FGF-2 on cell proliferation were observed.Neufeld and Gospodarowicz began with plates each containing 4 x 10 4 cells from a baby hamster kidney cell line (BHK-21).One set of plates was exposed to increasing concentrations of FGF-1 (ranging from 50 pg/mL to 250 ng/mL), while another set of plates was exposed to increasing concentrations of FGF-2 (ranging from 2.5 pg/mL to 25 ng/mL).These increments of growth factor were added in two boluses, one on day 0 and one on day 2.After 4 days, the number of cells on each plate was counted and recorded.These data were displayed as FIG 2 in [9], recreated here as Figure 1.
Experiment 2. Next, Neufeld and Gospodarowicz examined the ability of FGF-1 and FGF-2 to displace iodinated forms of both FGFs.This experiment functioned as a test to determine if FGF-1 could compete with FGF-2 the same cell surface receptor.In these experiments, the growth factors were labeled with 125 I. Then BHK-21 cell membranes were incubated for 30 minutes in the presence of either 125 I-FGF-1 (13 ng/mL) or 125 I-FGF-2 (14 ng/mL) and in the absence of competing unlabeled FGF.These values represented the maximal amount of 125 I-FGF-1 and 125 I-FGF-2 binding, respectively.Next, increasing concentrations of unlabelled FGF-1 and FGF-2 were added to 13 ng/mL 125 I-FGF-1 and then separately added to 14 ng/mL 125 I-FGF-2.Iodinated FGF binding was recorded and appeared in parts B and A, respectively, in figure 4 from [9], recreated here as Figure 2.2 Figures 1 and 2 serve as the standard against which we reference the fit of the models during development.We now begin formulation of the models by examining the biochemical kinetics involved in the competitive pathways of each experiment.

Experiment 1.
The first competitive pathway is as follows: Suppose R is a free receptor on a BHK-21 cell capable of being activated by either FGF-1 or FGF-2.Let G 1 be a molecule of FGF-1.Then, the binding of FGF-1 to a free receptor leads to an intermediate complex, {RG 1 }, which releases a product, call it P 1 , by the mechanism: (1)  The product P 1 begins a tyrosine-kinase signal transduction pathway leading to an increase in cell number. 3Concurrently occurring is the binding of FGF-2, G 2 , to another free receptor, R.This binding also leads to an intermediate complex, {RG 2 }, which again releases a product, P 2 : ( Product P 2 also initiates a signal transduction pathway.This cascade again results in increased cell numbers.When both species of growth factor are present, the interplay of these two equations, ( 1) and ( 2), results in competition of both species for the same free receptors: This competitive pathway is further developed by writing down the laws of mass action4 for (1) and (2), as follows 5 6 : At this point we employ the Michaelis-Menten steady-state assumption explained in [13].Essentially, this assumption states that following the initial stage of the reaction, termed the transient phase, the rate of synthesis of an intermediate remains approximately equal to the rate of consumption of said intermediate until the substrate, or growth factor in the present example, is nearly exhausted.Thus, a quasiequilibrium is reached.Applying this hypothesis, we take the concentrations of both intermediates to be constant, and using the notation where σ G 1 and σ G 2 are constants for cellular expression of G 1 and G 2 , respectively, and µ G 1 and µ G 2 are decay rates for the aforementioned growth factors, and [R] T is the total concentration of receptors.We may neglect σ G 1 and σ G 2 because the expression of either growth factor by the BHK-21 cells is negligible relative to the concentration of growth factor being added into the cell cultures.Likewise, we may neglect µ G 1 and µ G 2 because the decay rate of either growth factor is negligible relative to the concentration of growth factor being consumed by the growing cell populations.Essentially, FGF-1 and FGF-2 are being consumed by the cells at a far faster rate than either half-life would allow for decay. 6It is important to remark that the first and third equations in (4) have been simplified.Taking into account cell expression and FGF turnover rate, these equations are more completely written as: where σ G 1 and σ G 2 are constants for cellular expression of G 1 and G 2 , respectively, and µ G 1 and µ G 2 are decay rates for the aforementioned growth factors, and [R] T is the total concentration of receptors.We may neglect σ G 1 and σ G 2 because the expression of either growth factor by the BHK-21 cells is negligible relative to the concentration of growth factor being added into the cell cultures.Likewise, we may neglect µ G 1 and µ G 2 because the decay rate of either growth factor is negligible relative to the concentration of growth factor being consumed by the growing cell populations.Essentially, FGF-1 and FGF-2 are being consumed by the cells at a far faster rate than either half-life would allow for decay.
for each Michaelis constant, the second and fourth equations in (4) become: Substituting the first and second equations of ( 5) into the first and third equations of (4), respectively, yields: Next, we relate cell density to receptor concentration.We assume, as noted in [2], that the number of BHK-21 cells per unit volume is proportional to the total number of receptors that can initiate a signal transduction pathway in response to a growth factor.Thus, we may write where [N] denotes the concentration of BHK-21 cells, [R] T denotes the total concentration of receptors, and κ is the proportionality constant.Substituting R 0 /N 0 for the proportionality constant κ, we may write: where N 0 is the carrying capacity of the BHK-21 cells and R 0 is the total number of receptors at carrying capacity.As [2] explains, we may take R 0 to be on the order of unity; thus, our relationship becomes: Furthermore, we may write the total concentration of receptors as follows: Substituting the first and second equations of ( 5) into (10) yields: Solving for free receptors, [R], gives Substitution of ( 8) into (12) yields: [R] = [N ] Finally, ( 13) can be substituted into the first and second equations of ( 6), as follows: Describing cell proliferation is slightly more complex but accomplished when several biological considerations are taken into account.First, we assume that cell proliferation is logistic as determined from the characteristic shape of Figure 1.Secondly, as noted in [2], it is reasonable to assume that BHK-21 cell mitosis depends on the concentrations of both growth factors and BHK-21 cell apoptosis is linear in cell density.These considerations allow us to write: where φ(G 1 , G 2 ) is the coefficient of the logistic term and µ is the decay rate of BHK-21 cells.The term φ(G 1 , G 2 ) is a measure of how the growth factors influence mitosis.In the present model, φ(G 1 , G 2 ) takes the form: As explained in [2], the underlying idea is that sufficient concentrations of either growth factor are necessary for the birth rate to exceed the death rate, but the effects of FGF-1 and FGF-2 on birth rate at saturation of either growth factor are limited to a maximum value of λ.Thus, the equation for cell proliferation becomes: Combining this equation with the equations in (14), we obtain our first model, a system of three partial differential equations: Next we explore the competition of unlabelled FGF with iodinated FGF.The first competitive pathway of Experiment 2 is much the same as the pathway for Experiment 1 in that a free receptor, R, on a BHK-21 cell is capable of being activated by either of two species; however, in this pathway, the species are FGF-1 and iodinated FGF-2.Suppose as earlier G 1 is a molecule of FGF-1, and let G 2 * be a molecule of iodinated FGF-2.Then, the binding of FGF-1 to a free receptor leads to an intermediate complex, {RG 1 }, by the equilibrium: Again, there is competition for the free receptors; however, this time the competing species is iodinated FGF-2.The binding of iodinated FGF-2 to a free receptor also leads to the formation of the intermediate complex, {RG 2 * }, by the equilibrium7 : In a similar fashion we may write chemical equations for each of the remaining three pathways.Table 2 summarizes the species present in these pathways.

Table 2. Notation for Additional Species in Kinetic Equations
Species Notation fibroblast growth factor, FGF-1 We may write the dissociation expression for chemical equation (20) in terms of the intermediate, as follows: where the dissociation constant K 2 D is given by K The astute reader will notice that equation (21) takes the same form as the equations in (5), where now we are dealing with dissociation constants rather than Michaelis constants.Nevertheless, we may make a substitution for free receptors similar to that of equation (12) in experiment 1 to obtain a model for the first pathway in experiment 2, as follows: Using similar methods we obtain the models for each pathway in experiment 2. A summary of these four pathways is shown below: It is important to remark that the design of experiment 2, as carried out in [9], does not depend on the initial concentration of iodinated FGF nor the value of [R] T , because the final concentration of iodinated FGF is plotted as a fraction of the maximum.To understand this observation more fully, we examine the parameter [R] T .We see from ( 22 where [A] i represents the concentration of species A at each data point i.Here we see that, not only is the plot of [{RG 2 * }] independent of the initial concentration of iodinated FGF, it is also independent of [R] T and the entire numerator, including [G 2 * ], for we may simplify (24) to yield Hence, we rewrite the equations in (23), taking into account these observations of normalization, as follows: 4. Simulations and Optimization.Now that we have constructed a model for each of the competitive pathways, we use MATLAB to simulate the experiments performed by Neufeld and Gospodarowicz in [9].
4.1.Experiment 1.For the pathway involving differential equations, we use the MATLAB solver ODE15s for simulations.First, we simulate the initial trial performed by Neufeld and Gospodarowicz in [9].In this trial, cell plates containing 4 x 10 4 BHK-21 cells were exposed to increasing concentrations of FGF-1 in the absence of FGF-2.The added amounts of FGF-1 are shown in the first column of Table 3.

Table 4. Numerical Values of Parameters Used in Simulations
Parameter Numerical Value (from [2]) We now use the ODE15s solver to find the concentration of FGF-1 at time 48 hours. 8To this concentration of FGF-1 we add the second bolus of growth factor, again expressed in column 1 of Table 3.Finally, we use the solver to determine the number of BHK-21 cells at time 72 hours.
Using a similar method, we simulate increasing concentrations of FGF-2.In this trial there is no FGF-1 present, (i.e.[G 1 ] = 0), and the model does not depend on the values of k 2 and K 1 M from the equations in (18).Instead, this model utilizes the parameters µ, λ, N 0 , k 4 , and K 2 M .Numerical values for these parameters were again supplied by [2] and are given in Table 4. Again, by solving the system of differential equations twice, employing the pulse of additional growth factor described earlier, we obtain an approximation of the biological data.The data for both trials are plotted along with the associated biological data from Figure 1 in Figure 3. 9We now employ optimization to extrapolate the numerical values of the parameters appearing in the model.In the first trial of this experiment, these parameters are µ, λ, N 0 , k 2 , and K 1 M .To find the values of these parameters, which give the closest fit to the actual biological data, we first define an error function.This function is the sum of the squares of the differences of the biological data for cell density and the data calculated from the model for cell density, as represented below: This error function has the values of the parameters as inputs.Different values for the parameters yield a different numerical value for the error function.Then, using a tool in MATLAB known as fminsearch, we minimize the error function and the resultant output is a vector of the values of the parameters which give the closest fit to actual biological data.Using fminsearch for the first trial, the resulting coefficient vector is: 0.014162 0.468937 790, 568 1.952589 0.012074 , which corresponds to the values for µ, λ, N 0 , k 2 , and K 1 M .Likewise, we apply an error function to the second trial.Here, we are searching for the values of the The revised model, taking into account the optimal values of the parameters, is plotted along with the accompanying biological data in Figure 4. Now we must consider the overlap between the two trials.The individual optimizations yielded slightly different values of µ, λ, and N 0 , as shown in Table 5.

Table 5. Comparison of Shared Parameters
Parameter FGF-1 trial FGF-2 trial µ 0.014162 h −1 0.019587 h −1 λ o.468937 h −1 0.775144 h −1 N 0 790,568 cells 789,977 cells Using a combined error function where the parameters are defined only once should give a compromise fit for the two trials.This combined error function yields the coefficient vector: 0.020241 0.64678 810, 667 1.3847 0.02925 0.077062 0.012556 , which corresponds to parameters µ, λ, N 0 , k 2 , K 1 M , k 4 , and K 2 M .Figure 5 shows a plot of the model utilizing these parameters.4.2.Experiment 2. Next we use MATLAB to simulate the second experiment performed by in [9].Recall that BHK-21 cell membranes were incubated for 30 minutes in the absence of competing unlabelled FGF and in the presence of either 13 ng/mL 125 I-FGF-1 for part A or 14 ng/mL 125 I-FGF-2 for part B. To this  concentration of iodinated FGF, various concentrations of unlabelled FGF-1 and unlabelled FGF-2 were added, as shown in Table 6.These added concentrations of growth factor along with values for the parameters K 1 D and K 2 D , which were approximated in [9] and are displayed in Table 7, allow us to calculate the concentration of iodinated FGF for each pathway in equations (26).A plot showing the model for each of the pathways in parts A and B appears along with the associated biological data in Figure 6.
We now employ optimization to extrapolate the numerical values of K 1 D and K 2 D .In order to find the optimal values of these two parameters, we again define an error function.This function is the sum of the squares of the differences of the ordinate from the biological data, p exp i , and the ordinate calculated from each pathway, p model i , as represented below: Using fminsearch for this error function, the resulting coefficient vector for the parameters K 1 D and K 2 D is: 3.734 x 10 −4 6.65 x 10 −5 .
The revised model for parts (A) and (B), taking into account the optimal values of K 1 D and K 2 D , is plotted in Figure 7 along with the associated biological data.
we may rewrite K 1 M as follows: Thus, k 1 is given by: Substituting in the values of k 2 , K 1 M , and K 1 D , yields a value for k 1 of 47.9507.Moreover, we may solve the equation of K 1 D for k −1 as follows: Substituting in the values of K 1 D and k 1 gives a value for k −1 of 0.0179.Likewise, we can compute the values of k 3 and k −3 .We find that k 3 is 6.1699 and k −3 is 0.00041.Table 8 displays these values along with the other optimized parameters.M and K 2 M as previously thought.This result has a very important implication.It shows that k 2 and k 4 are not always insignificant, and this fact must be taken into consideration before simply disregarding the values of these parameters.Moreover, the fact that k 2 and k 4 are significant greatly affects the difference between K i M and K i D .Our results show that K 1 M and K 2 M are fifty to one hundred times greater than K 1 D and K 2 D , respectively 10 .This result again has implications for future research. 10It is important to note that the values of K 1 D and K 2 D predicted by our model, shown in Table 8, show good agreement with the values estimated in [9] and shown in Table 7.
Second, our model demonstrates a close relationship between FGF-1 and FGF-2.By taking the ratio between the first and second equations of ( 14), we may write Separating the variables, integrating, and exponentiating equation (34), yields Thus, we have shown that either growth factor can be represented in terms of the other, if the initial values of both are known.This explains why the light blue and dark blue curves in the Figures 6 and 7 are parallel.Finally, our model demonstrates the importance of parameter estimation in the modeling of biological phenomena.Slightly changing the values of parameters embedded in mathematical models can result in noticeable changes in the fit of the model.This result was shown with the optimizations we performed.Now that we have constructed our model and the parameters have been accurately estimated, we can use the model to make predictions.From our model we are able to formulate several testable hypotheses using MATLAB.We hypothesize about the appearance of several variations of experiment 1.First, we predict the outcome if the number of pulses is changed.Figure 8 shows the results of a replication of experiment 1, the only difference being that in the first trial growth factor is added only initially, the second trial follows the exact procedures of experiment 1, and the third trial entails a pulse of growth factor added initially and consecutively each of the next three days.Each of these trials still involved counting cell number at four days, and the same total amount of growth factor was added for all three.Only the average pulse size was varied for each trial.
Next, we predict the outcome of changing the number of days we wait before counting.Figure 9 compares the start of counting on day 4 instead day 6, when pulses are added initially, day 1, day 2, and day 3.
Here we notice that the greatest cell density occurs when counting earlier (day 4) rather than waiting to count (day 6).This could be attributable to the decay of the growth factor of the BHK-21 cells.
Next, we compare the effects of adding growth factor on consecutive or alternating days and then counting on day 6, as shown in Figure 10.
Here we observe that the model predicts that the trials will initially overlap and then diverge later in the experiment.This could be attributable to the growth factor decaying in the alternating trial, while the consecutive trial has enough growth factor to last longer.
Finally, we predict the effects of adding growth factor for an increasing number of pulses.As Figure 11 shows, successive trials attain greater cell density when growth factor is added for a greater number of days.Cells/dish FGF−1; day 0 FGF−2; day 0 FGF−1; day 0, 2 FGF−2; day 0, 2 FGF−1; day 0, 1, 2, 3 FGF−2; day 0, 1, 2, 3  Here we have demonstrated that the interaction of multiple growth factors with cell surface receptors can be modeled to produce predictable outcomes.Our model correctly describes the results of experiments performed in [9] and can predict the outcome of many experimental protocols, given accurate parameters for modelling.Although the current model was developed to simulate a relatively simple cell culture system with only two growth factors and one receptor, its capacity for expansion to include more growth factors and growth factor receptors identifies this model as an excellent base for developing testable simulations of complex biological systems.The development of predictive models is essential to understanding the complex interplay of growth factors and their receptors, as happens during embryonic development and wound healing.

Figure 2 .
Figure 2. Competitive Inhibition in the Binding of 125 I-FGF-1 and 125 I-FGF-2 to BHK-21 Cell Membranes by Unlabelled FGF-1 and Unlabelled FGF-2 (from [9]) ) that the concentration of 125 I-FGF-2 depends on [G 1 ], [G 2 * ], K 1 D , and K 2 D .Furthermore, we note that only the value of [G 1 ] changes and thus [{RG 2 * }] varies only with changes in [G 1 ].Now we consider the normalization of [{RG 2 * }].Each entry of [{RG 2 * }] is divided by the first entry, the maximum, to yield a fraction of the total binding.Thus, we may write (22) as follows:

2 Figure 3 .
Figure 3. Initial Fit of Model to Biological Data from [9].

Table 1 summarizes the species present in this pathway.Table 1 .
Notation for Species in Kinetic Equations

Table 6 .
Concentrations of FGF-1 and FGF-2 Added in Parts A and B of Experiment 2

Table 7 .
Numerical Values of Additional Parameters Used in Simulations

Table 8 .
Numerical Values Extracted from Optimizations Discussion and Future Work.A number of findings can be drawn from our model.First, our model gives a new perspective on the relationship between K i M and K i D .Our model demonstrates that k 2 and k 4 are the driving force for experiment 1, and not K 1