A mathematical model quantifies proliferation and motility effects of TGF--$\beta$ on cancer cells

Transforming growth factor (TGF) $\beta$ is known to have properties of both a tumor suppressor and a tumor promoter. While it inhibits cell proliferation, it also increases cell motility and decreases cell--cell adhesion. Coupling mathematical modeling and experiments, we investigate the growth and motility of oncogene--expressing human mammary epithelial cells under exposure to TGF--$\beta$. We use a version of the well--known Fisher--Kolmogorov equation, and prescribe a procedure for its parametrization. We quantify the simultaneous effects of TGF--$\beta$ to increase the tendency of individual cells and cell clusters to move randomly and to decrease overall population growth. We demonstrate that in experiments with TGF--$\beta$ treated cells \textit{in vitro}, TGF--$\beta$ increases cell motility by a factor of 2 and decreases cell proliferation by a factor of 1/2 in comparison with untreated cells.


Introduction to the Biology of TGF-β
In normal organisms, the growth of cells is under tight regulation by growth factors and is highly dependent on the developmental stage during the lifespan of the organism. Disruption of this regulation is the most frequent cause of cancer diseases. Unlike normal differentiated cells, cancer cells are usually hyperproliferative as a result of the abnormal activation of multiple growth-stimulating intracellular signalling pathways and loss of tumour suppressors. In cancer cells, these pathways do not respond to normal regulatory signals but are manipulated by one or more oncogenic signals, often encoded by oncogenes. Expression of oncogenes alters signalling pathways that under normal conditions maintain cell growth homeostasis. Thus, altering these pathways may favour increased cell and cancer growth. A good example is the transforming growth factor (TGF) β family, which is known to be able to act as both a tumour suppressor and tumour promoting factor.
The TGF-β family consists of multitasking cytokines that play important roles in cell proliferation, cell motility, apoptosis, lineage determination, extracellular matrix production, and modulation of immune function [22]. These ligands bind to a heteromeric complex of transmembrane serine/threonine kinases, the type I and type II receptors (TβRI and TβRII). The receptors are activated upon ligand binding, leading to the subsequent phosphorylation and activation of a family of transcription factors called Smads, which regulate transcription of a subset of genes [23]. In addition to Smads, other signalling pathways have been implicated in TGF-β actions in recent studies. These include the extracellular signal-regulated kinase (ERK, MAPK), c-Jun NH2-terminal kinase (JNK), p38MAPK, phosphatidylinositol-3 kinase (PI3K), and Rho GTPases (reviewed in [7,11,39,43]). The critical role of these non-Smad pathways on mediating the cellular effects of TGF-β remains to be fully characterised.
TGF-β was originally reported to induce transformation of mouse fibroblasts [24]. Subsequent studies indicated that TGF-β is a potent inhibitor of cell proliferation and a tumour suppressor [32,36]. Consistent with its tumour suppressor role, many cancers lose or attenuate TGF-β-mediated anti-mitogenic action by mutational inactivation of TGF-β receptors or their signal transducer Smads [13,15,16,20,40,41]. There is increasing evidence to show that excess production and/or activation of TGF-β in tumours can accelerate cancer progression through enhancement of tumour cell motility and survival, increase in tumour angiogenesis, extracellular matrix production and peritumoural proteases, and the inhibition of immune surveillance mechanisms in the cancer host (reviewed in [7,9,11]). Cancer progression and metastasis consist of a series of sequential events. After initial cell transformation, often mediated by the function of oncogenes, tumour cells growing at the primary site will invade the surrounding stroma and migrate towards blood vessels. Through various mechanisms such as epithelial-mesenchymal transition (EMT), tumour cells will enter the blood vessels and travel to other parts of the body through the circulatory system. Some of the cells will then arrest at distant sites where they may proliferate and invade into the adjacent organs. Cell motility is therefore a critical element during the spread of tumour cells from their initial sites of residence. In this study, we focus on the tumour-promoting effect of TGF-β through inducing cell motility.
The receptor tyrosine kinase HER2 (ErbB2, Neu) belongs to the family of epidermal growth factor receptor (EGFR). Gene amplification or overexpression of HER2 is observed in about 25% of breast cancers. TGF-β has been shown to synergize with the oncogene ErbB2 in cancer progression. Overexpression of active TGF-β 1 or active mutants TβRI (Alk5) in the mammary gland of bigenic mice also expressing mouse mammary tumour virus (MMTV)/Neu (ErbB2) accelerates metastases from Neu-induced mammary cancers [25,26,27,35]. Exogenous as well as transduced TGF-β confer motility and invasiveness to MCF10A nontransformed human mammary epithelial cells (HMEC) stably expressing transfected HER2 [34,38]. Expression of the oncogene HER2 in these cells does not affect the function of TGF-β on inhibiting cell proliferation [38]. It is likely that in many cancers, TGF-β may still attenuate proliferation while inducing cellular events associated with metastatic dissemination, such as cell motility. In this paper we report experiments with MCF10A/HER2 cells to study and to separate the effects of TGF-β on cell proliferation and motility. Due to the complexity of TGF-β signalling that simultaneously affects several biological parameters, it is important to computationally simulate the behaviour of cells under TGF-β exposure. Our model, which can also be adopted to simulate other growth-regulating signals, will provide a unique insight into the TGF-β function in both normal and cancer cells, and further understanding of targeted therapeutic strategies that aim at interfering with TGF-β signalling.
Mathematical modeling of chemotherapy with Tamoxifen has been carried out in [37], with special consideration of TGF-β. There, the authors use a discrete cellular automaton (CA) model to investigate the effects of TGF-β on tumour morphology and invasiveness. It is possible to assign properties to individual cells in a cellular automaton model (such as different levels of malignancy). However, they can be computationally costly, in particular in situations where large cell populations are involved. Since we are interested chiefly in the effects of TGF-β on a large number of cells growing in a plate, rather than individual cell behaviour, we work with a continuous partial differential equation model. Future research problems may require the construction of hybrid models such as in [1] where discrete cells are coupled to continuous fields that in turn are solutions of partial differential equations.

Introduction to the Mathematical Model
The Fisher-Kolmogorov equation (or KPP equation) has been studied and used widely in mathematical biology. It originated in the 1930s in works of R. A. Fisher and A. N. Kolmogorov [12,18] where it was proposed as a model for the spread of an advantageous gene in a population. Mathematically, it falls into the large class of reaction-diffusion equations, in which one or more diffusible "species" enter into a scheme of reactions. The species can be chemical substances or biological populations. Reaction-diffusion equations in general and the Fisher-Kolmogorov equation in particular have been applied successfully by many authors to model the behaviour of cells in tissues, see for example [1,2,5,6,8,10,19,30]. The Fisher-Kolmogorov equation is usually used to show traveling wave behaviour, but we use it here to model the spatial growth of a cell culture from initial seeding to confluence. Let u(x, y, t) denote the spatial density of tumour cells at time t in a dish corresponding to a spatial region Ω ⊂ R 2 . The time evolution of this density is described by the partial differential equation ∂ ∂t u(x, y, t) = D∆u(x, y, t) + αu(x, y, t)(1 − βu(x, y, t)).
Here α incorporates the intrinsic proliferation rate (unit time −1 ), from initial seeding to confluence, in a single value independent of x, y, t. Also, β denotes the area occupied by an average single cell. The density satisfies, u(x, y, t) ≤ β −1 for all x, y, t, and the total number of cells in the region Ω ⊂ R 2 at time t is given by The spatial region Ω may be the entire region occupied by the population or a subfield thereof. Equation (1) is accompanied by an initial condition and appropriate boundary conditions. If a subfield Ω 1 of Ω is completely populated, then u(x, y, t) = β −1 for all x, y ∈ Ω 1 , otherwise u(x, y, t) ≤ β −1 and u(x, y, t) dx dy is the fraction of possible cells in Ω 1 at time t. Here and in the sequel | · | denotes the area of the domain in question. We denote the total mass of the density by (2) Much of the remainder of this paper will be concerned with the problem of finding parameter values to a modified version of equation (1) from experimental observations and with the interpretation of the solution to the equation. The diffusion constant D in equation (1) has to account for three simultaneous processes contributing to the spatial movement of cells. First, dividing cells occupy increasing space, therefore a certain part of D, say D p , has to account for the spatial expansion of proliferating cells. The authors of [6,8,19,30] have proposed relations of the type Here is related to the typical length or diameter of a cell and T d is the doubling time. The dimensionless factor k 1 in [8] accounts for the expansion of the viable rim of dividing cells. We take to be the average increase in cell diameter between mitotic events. We assume the area of an average cell increases from approximately 50 µm 2 to 100 µm 2 , so we take = 3.3 µm. We assume the average cycle time is . Second, individual cells plated on a Petri dish do not remain fixed in their position but undergo a random walk. It is natural to assume that the random walk is not biased and follows the laws of Brownian motion. Third, cells in a cluster may break lose from that cluster as TGF-β is known to decrease cell-cell adhesion [28]. Thus, a part of D, say D m , must account for this second and third contribution to cell movement. We propose here a modification of the standard Fisher-Kolmogorov equation (1). Our reasoning is that as the cells become more densely packed, their random motility should decrease. Thus, with U from equation (2) we propose The constant D m , which captures the random motility of individual cells and clusters and the tendency of clusters to break apart, depends on TGF-β concentration.
We remark here that D could also be made dependent on the local density u. In the case of interest to us in the present paper the initial seeding of the cells results in a uniform dispersion of the population that justifies the approximation made in equation ( Table 1-1] D m = 5 · 10 −10 cm 2 s −1 = 180 µm 2 h −1 . Chaplain and Matzavinos [6] arrive at a similar value D m = 7 · 10 −5 cm 2 day −1 = 300 µm 2 h −1 . This value is obtained in [6] from the Einstein- where k B is Boltzmann's constant, T the temperature, d the diameter of the cells, and η the dynamic viscosity of the surrounding medium. It can be excluded, however, that cells in vitro, such as in our experiments, move in an environment that has the viscosity of water, η H2O = 10 −3 P a s. Thermally driven motion and Stokes friction are not relevant in the context of moving cells.

Materials and Methods
MCF10A/HER2 cells were generated and grown as described [38]. All live cell imaging experiments were performed in triplicate. For growth assay, equal number of cells (1.5 · 10 4 /well) suspended in full medium were seeded on 6-well plates. Cells were allowed to grow in the absence or presence of TGF-β at different concentrations (0.5 − 5 ng ml −1 ) over a time course of 8 days. Medium containing fresh TGF-β was replenished every 2 days. Cells were harvested by trypsinization every 24 h, and subjected to total cell number counting using a Coulter counter.
In Figure 4 B each data point represents the mean of 3 wells. Standard deviations were less than 5% in all cases and were not included in the figure. For random motility assay of individual cells, cells were seeded on 6-well plates at 3 · 10 4 /well one day prior to image recording. For motility of cell clusters, cells were seeded at 10 4 /well two days prior to image recording. Culture media was replaced with fresh serum-free or full L-15 medium and different concentrations of TGF-β was added 1 h before the image recording. Cells were imaged on a Nikon TE 2000-E microscope and captured using a Hamamatsu Orca ER camera and Metamorph software. Images were taken every 5 minutes for 270 minutes for random motility assays and every 30 minutes for 14 hours for the cell cluster experiments. Image analysis was performed using Metamorph software (MDS, Inc., Toronto, Canada).

Random motility of individual cells
We have taken position measurements of individual cells and have obtained time series data (x i , y i ) N i=1 , with ∆T = 5 min. The mean-squared displacement (MSD) is calculated according to where k = 1, 2, . . . , k max . The maximal multiple k max for which r 2 is calculated has to be chosen carefully. Obviously, the bigger k, the fewer such displacements are contained in the time series of length N thus the larger the uncertainty in the estimate r 2 (k). The paper [31] gives estimates of the standard deviation of r 2 (k) in terms of N and k max . To obtain the random motility of a single cell we fit the estimate (4) with a linear function (see [3] for a discussion of this relation) Unfortunately, the analysis is made difficult by cells that seem to remain immobile or whose r 2 (k)-curves do not look linear. For the plots in Figure 1 we have selected only curves that could be fit reasonably by a straight line. Other curves show cells moving initially and then coming to a standstill. Moreover, we cannot see a clear increase of D m as TGF-β is added. For the average random motility of individual isolated mobile cells we obtain, both for cells treated with 0 and with 5 ng ml −1 TGF-β, This value can be taken as an upper bound which is realised only if all cells are mobile.
We now look at the percentage of mobile cells. We define a cell to be mobile if it moves outside of a 100 µm × 100 µm square centred at the cell's original position. It becomes clear ( Table 1) that the percentage of mobile cells in our experiment increases as TGF-β is introduced.
It seems to be a frequently observed phenomenon that a certain subpopulation of cells remains immobile. Selmeczi et al. [33] observed different types of keratinocytes and fibroblasts and selected trajectories only for cells that were moving at all. This selection procedure is clearly admissible for single particle tracking experiments, but is not applicable to a population with a high fraction of non-mobile cells. We want to emphasize that the cells are not subject to a directional gradient of TGFβ, rather a uniform concentration applied from all directions, therefore they are able to move randomly in any direction, as opposed to being stimulated to move in one particular direction. Although it is unknown why cells suddenly stop moving, it is possible that in the absence of a strong directional signal there is "noise" in the motility signaling that leads to immotility in randomly moving cells.

Random motility of cell clusters
We next investigated the motility of cell clusters (≈ 10 cells) in the absence or presence of TGF-β. As shown in Figure 2, top row, untreated cell clusters are relatively immobile and seem to maintain a certain shape, although the cells can still rotate within the clusters. In contrast, clusters treated with TGF-β are significantly more active in altering their shape and migrating towards unoccupied space ( Figure 2, bottom row). We analysed the movie data in the following way. On the initial frame of the movie the centre of the cluster is identified as the centre of the smallest rectangle containing the cluster. Then, on each subsequent frame of the movie, a straight line is drawn to the point on the boundary of the cluster the farthest from the initial centre. The time gap between two frames is ∆T = 1 h. We obtain a star of radii as shown in Figure 3 A. For a time series of N radii (r i ) N i=1 we calculate the variation of the squares For a total of 36 movies, 18 control cases and 18 treatment cases (5 or 10 ng ml −1 ) respectively, a difference can be seen in the value of the variation v 2 (Figure 3 B) with TGF-β treated cells showing a higher value than untreated cells. Movies were obtained both under serum-free conditions and with serum. The clusters treated with TGF-β are more active by a factor of two.
We combined the observations of individual mobile cells and clusters of cells to obtain values of D m = 5 µm 2 h −1 (TGF-β concentration 0) and D m = 10 µm 2 h −1 (TGF-β concentration 5 ng ml −1 ). These values are scaled downwards from the upper bound obtained from the single cell motility experiments in view of the fact that only a fraction of cells are fully mobile. The ratio of the D m 's for treated and untreated case is suggested by our experimental observations.

Growth assays and their numerical simulation
The mathematical model from equation (3) was programmed using matlab (version 7.1, The MathWorks, Inc., Natick, MA). A standard Crank-Nicolson scheme is used for the diffusion step [29,Section 17.3], combined with an Euler forward reaction step. The codes will be available upon request from the corresponding author.
We fix the following length and time scales for our simulations. The square domain has side length L = 300 µm and the time unit is T = 1 h. We assume that a single cell occupies an area of β = 100 µm 2 , thus the field has a carrying capacity of 900 cells. We impose periodic boundary conditions to account for the fact that our numerical domain is a small section of a larger domain, say a Petri dish. Thus, a cell mass that leaves the region will be balanced by a cell coming in at the opposite side. It should also be pointed out that the carrying capacity is the same in control and treatment cases.
Numerical simulations of cells growing in a field are shown in Figure 4 A. The initial datum resembles the random placement of cells on the Petri dish. Using our modified Fisher-Kolmogorov equation (3) we can simulate the experimentally determined growth data (Figure 4 B, discrete symbol represent experimental data, continuous curves represent simulations). At first, the experimental growth data were normalised with respect to a maximum capacity of 1.1 · 10 5 cells. The parameter α was obtained from fitting the experimental data (see Fig. 4 B), while D m and D p were chosen as stated above. The numerical simulation starts only 72 h into the experiment. In our previous experiments with these cells [17], we have noticed that after the cells suspended in medium are seeded onto the plates, they require about 48 hours to adhere to and spread on the plates and adjust to their new environment, before they start a typical growth regime. The numerical parameter values are collected in table 2. To support our choice of α we note that if the cells were able to grow exponentially according to the law u = αu, then the corresponding doubling times T 2 = ln 2/α would be T − 2 = 10 h for untreated cells (α = 0.07 h −1 ) and T + 2 = 17 h for cells treated with TGF-β (α = 0.04 h −1 ). These doubling times are feasible for MCF10A/HER2 cells in the early exponential growth phase.

Discussion and Conclusions
Our experiments demonstrate that TGF-β increases the percentage of mobile cells in a cell population in a dose-dependent manner, rather than increases the mean square migration displacement of individual cells. As we demonstrated, random migration of cell clusters containing ≈ 10 cells is a feasible approach to parametrise the unbiased random cell migration in a large population of cells. The cluster motility assay indicates that clusters are more mobile and/or less cohesive if TGFβ is present. In fact, biological evidence for both mechanisms has been identified. TGF-β induces a PI3K-mediated activation of Rac1/Pak1 pathway, which leads to increased cytoskeleton reorganisation, turnover of focal adhesion, and eventu-ally cell migration [42]. Furthermore, TGF-β downregulates cell-cell adhesion by decreasing the adhesion protein E-cadherin [28]. The cluster motility assay reflects an overall effect of the multifunctional molecule TGF-β on motility. In metastatic cancer cells, TGF-β stimulates epithelial-mesenchymal transition (EMT) which is necessary for the acquisition of invasive and metastatic phenotype. It has been shown that in cell monolayer wounded with a pipette tip, the presence of TGF-β in the medium markedly induces the closure of wounded area, which is another indicator of TGF-β-induced cell motility. However, this assay is not used in our parametrisation as the cell migration is not random but targeted to the wounded area.
Our model has a relatively small number of parameters: D p , β, D m , and α. Our procedure for determining these parameters is as follows: (1) D p and β are determined by the average doubling time and average area of cells in the culture, respectively. (2) D m is determined by adjusting the motility measurements of isolated mobile cells and clusters of mobile cells for their fraction of the total population as the culture attains confluence. From our experiments the upper estimate on D m was 120 µm 2 h −1 , which we adjusted downward by ≈ 1 order of magnitude to compensate for the fraction of mobile cells in the total population. Table 1  It is possible that different values for D m and α give rise to very similar total population growth curves. We determine ratios of D m and α for treatment and control cases that simulate our experimental growth data (Figure 4 B). Our determination of D m (corresponding to motility) and of α (corresponding to proliferation) reveals that TGF-β roughly doubles the value of D m and halves the value of α. The total population growth data of the cell cultures alone (without any further information) are insufficient to quantify these two separate effects of TGF-β. We have used in addition the motility experiment data to distinguish the two effects.
We have developed a general spatial model of proliferating cell cultures in vitro, which allows quantification of the properties of proliferative capacity, cell mobility, and clustering as the population attains confluence. The novelty of our interpretation of the Fisher-Kolmogorov equation is that cells begin as isolated geometric regions corresponding to seeding. As they proliferate and migrate these regions expand as controlled by the diffusion parameters and proliferation parameter. When the density is identically β −1 (whose value is readily correlated to the average area occupied by an individual cell), the region occupied is completely populated. When the density is < β −1 in a region, then the integral of the density over this region is the expected value of the number of cells in the region. Our focus here is on quantification of the effects of TGF-β on motility and proliferation for in vitro cancer cell cultures, but our model allows investigation of similar phenomena for many cell types, both prokaryotic and eukaryotic. Extensions and modifications of our model could incorporate different cell proliferation dynamics, multiple cell types and cell interactions (for example endothelial cells, hematopoietic cells, and immune cells), in both in vitro and in vivo settings. In the future we plan to add more components that naturally exist in the tumour microenvironment to the model, including stromal fibroblasts, immune cells and vascular endothelial cells, which are all affected by TGF-β. In addition, many tumours carry key mutations on oncogenes or tumour suppressor genes, which may affect their responses to TGF-β. For example, cancer cells can become unresponsive to the anti-proliferative function of TGFβ, while remaining sensitive to TGF-β-induced motility. Inactivating mutations of TGF-β receptors are frequently found in pancreatic, biliary and colon cancers [14,21,40]. Usually, only a subset of cancer cells in a tumour harbour these mutations and these cells live together and communicate with other wild-type epithelial cells, fibroblasts and endothelial cells. It will be important to further investigate with the proper assistance from a mathematical model how these cancers are affected by TGF-β. The dual function of TGF-β as both a tumour suppressor and promoting factor is reflected by its individual effects on inhibiting cell proliferation and stimulating cell motility. Primary tumour growth by active cell proliferation results in the increase of tumour size. However, during cancer progression, dissemination of tumour cells from primary site into circulation and the subsequent invasion of other organs, in which cell motility plays a critical role, usually determine the grade of malignancy and the outcome of cancer treatment. Therefore, despite the fact that TGF-β has been considered a promising therapeutic target to treat certain types of cancer, the uncertainty of the overall effects of TGF-β intervention due to the potential risk of derepressed tumour cell proliferation has been an unsolved critical issue. This requires further in-depth understanding of the different aspects of TGF-β action through highly quantitative approaches such as mathematical modeling. Our model described here is the first step of encapsulating the multiple functions of a signalling molecule and potential therapeutic target such as TGF-β into a computational model. By introducing the modeling method to the cancer biology of TGF-β, questions that have been raised for years will have the chance to be answered from a novel angle. These interesting topics include: (1) When is the critical time of TGF-β being switched from a tumour suppressor to a promoter? (2) When is the most efficient time to apply anti-TGF-β therapy during cancer progression? (3) Will therapeutic inhibition of TGF-β at pre-cancer lesions or at early stages of the disease trigger primary tumour growth as a result of derepression of cell proliferation? Besides the complexity of TGF-β signalling, cancer cells can respond differently to TGF-β. For example, several oncogenes including HER2 are known to synergize with TGF-β to function on cell motility and survival. Some cancers do not respond to the anti-proliferative effect of TGF-β but still respond to its stimulatory effect on motility. The mathematical model will therefore provide a unique way to simulate how TGF-β functions differently in cancer cell populations with different characteristics, and the effect of TGF-β therapeutic inhibition on a tumour consisting of heterogeneous cell populations.      -TGF-β + TGF-β Figure 3. A Coverage radii for cell clusters in absence (blue) and presence (red) of TGF-β. The fact that all vectors point into the first quadrant is due to a minor technicality. B The variation of the squared radii is computed according to equation (6) for six different movies. The control cases are shown in blue, the treatment cases in red. These figures are representative for three independent realisations of the experiment.