A Macroscopic Mathematical Model for Cell Migration Assays Using a Real-Time Cell Analysis.

Experiments of cell migration and chemotaxis assays have been classically performed in the so-called Boyden Chambers. A recent technology, xCELLigence Real Time Cell Analysis, is now allowing to monitor the cell migration in real time. This technology measures impedance changes caused by the gradual increase of electrode surface occupation by cells during the course of time and provide a Cell Index which is proportional to cellular morphology, spreading, ruffling and adhesion quality as well as cell number. In this paper we propose a macroscopic mathematical model, based on advection-reaction-diffusion partial differential equations, describing the cell migration assay using the real-time technology. We carried out numerical simulations to compare simulated model dynamics with data of observed biological experiments on three different cell lines and in two experimental settings: absence of chemotactic signals (basal migration) and presence of a chemoattractant. Overall we conclude that our minimal mathematical model is able to describe the phenomenon in the real time scale and numerical results show a good agreement with the experimental evidences.

Introduction

Despite significant progress regarding potential therapeutic targets aimed at improving survival, patients affected by solid tumours frequently die for systemic spread of the disease to distant sides. Indeed, when cancer cells acquire the ability to separate and move away from the primary tumour mass, migrate through the surrounding tissue, and enter the lymphatic system and/or blood circulation, the prognosis becomes poor. Therefore, the control of cell motility is a new and attractive approach for the clinical management of metastatic patients. The quantitative assessment of tumour cell migration ability for each patient could provide a new potential parameter predictive of patient outcomes in the future.

To metastasise, tumour cells have to early acquire the ability to move and respond to motogen gradients [1]. Cell migration is a spatially and temporally coordinated multistep process that orchestrates physiological processes such as embryonic morphogenesis, tissue repair and regeneration, and immune-cell trafficking [2]. When cell migration is deregulated, it contributes to numerous disorders including tumour metastasis [3, 4]. Due to its important role in regulating physiological and pathological events, methods aimed to examine cell migration may be very useful and important for a wide range of biomedical research such as cancer biology, immunology, vascular biology, and developmental biology. Migrating cells respond to a plethora of mitogen stimuli, and serum (as mixture of growth factors, cytokines and chemokines) is a major source of chemoattractants. These chemoattractants, through the interaction with their cognate receptors allow cells to acquire a polarized morphology with the extension of adhesive protrusions [4]. This is followed by the attachment of the protrusion to the substratum at the cell front, the translocation of the cell body and, finally, the detachment of the trailing end of the cell from the substratum [5, 6]. Such a complex process requires the coupling of extracellular signals with the internal signalling machinery that controls cytoskeleton dynamics [7].

The most widely used technique to study cell motility in vitro is the Boyden chamber assay in which cells placed in the upper compartment of the chamber are allowed to migrate through a microporous membrane into the lower compartment, in which chemotactic agents are present; after an appropriate incubation time, the membrane between the two compartments is fixed, stained, and the number of cells that have migrated to the lower side of the membrane is determined [8]. The subjective nature of measurements and the inability to assess cell motility along the time are the major limitations of this assay.

Current molecular studies are providing a more global physicochemical picture of cell locomotion in which the role of spatial and temporal components of the process are detailed [9]. Recently, to overcome the manual and highly subjective nature of measurements, accelerate analysis and translate conventional Boyden chamber assay into an automated, quantitative high-throughput system, ACEA Biosciences developed the xCELLigence Real Time Cell Analysis (RTCA) technology able to automatically monitor cell motility in real-time without the incorporation of labels. The xCELLigence RTCA technology measures impedance changes in a meshwork of interdigitated gold microelectrodes located at the bottom side of a microporous membrane (CIM-plate). These changes are caused by the gradual increase of electrode surface occupation by migrating cells during the course of time and provide an index of cell migration. The relative electrical changes during a measurement are displayed by xCELLigence software as a unit less parameter termed Cell Index, which is calculated as a relative change in actual impedance divided by a previously registered background value. This method of quantitation is directly proportional to cellular morphology, spreading, ruffling and adhesion quality as well as cell number [10, 11]. To reach a quantitative understanding of the mechanisms underlying these processes, concepts and methods from mathematics and physics can be extremely valuable, as we will see in the following.

In general, mathematical models can be very useful to modelize a wide variety of biological systems including cell dynamics and cancer [12-18]. In particular, the development of quantitative predictive models, based on biological evidence, whose parameters are calibrated on biological data, can help in saving time and resources when designing novel experiments. Moreover, even though a mathematical model is not aimed to replace a real experiment, it can represent a guide to interpret acquired biological data and investigate new insights. In relation to in vitro assays in tumour chamber some mathematical model have been already proposed in the scientific literature, mainly focused on cell invasion experiments. In such context, the invasive ability of the cells is measured by the placement of a coating of extra-cellular matrix proteins on top of the porous membrane. In the papers [19, 20] a continuous model, based on partial differential equations (PDEs), was proposed in relation to a Boyden like cell invasion experiment. Then, in [21] the authors proposed a similar model to investigate some modulating factors of the cancer cell invasion, making also use of a real-time impedance-based in vitro technology.

In this paper, first we analysed basal (absence of chemotactic gradient) and directional (presence of serum as a source of chemotactic agents) cell migration of three different cell lines by the xCELLigence cell analyser: Melanoma A375, fibrosarcoma HT1080, and chondrosarcoma Sarc cell lines have been previously characterized for their migration ability by us and used as models of three different tumour types [22-24]. Then, we apply a theoretical analysis to describe cellular motility events, and we propose a mathematical model for the cell migration assay using the xCELLigence Real Time Cell Analysis (RTCA) technology. The proposed macroscopic model, based on advection-reaction-diffusion equations, adapts and extends the mathematical models in the aforesaid cited papers to the specific in vitro experiment in our analysis (see section Discussion and conclusions). We calibrated model parameters using real data, as well as information available in scientific and modellistic literature. With such estimate we carried out numerical simulations to compare simulated behaviour with the experimental data in absence or presence of chemotactic gradient. Our numerical results show a very good concordance with the experimental curves. Finally, we validated the model, simulating different experimental conditions, as the initial cell density, and then comparing numerical curves with data obtained from relative experiments. In this regard recorded experimental data on chondrosarcoma Sarc cell line seem to confirm theoretical results.

Results

Basal and directional cell migration of three different cell lines

For this study we considered three human, neoplastic cell lines which we have previously characterized for their cell migration ability [22-24]. Melanoma A375, fibrosarcoma HT1080, and chondrosarcoma Sarc cell lines (Fig 1 panel A) were subjected to both cell proliferation and migration assays using the xCELLigence Real Time Cell Analysis (RTCA) technology. This technology measures impedance changes in a meshwork of interdigitated gold microelectrodes located at the well bottom (E-plate for proliferation assay) or at the bottom side of a microporous membrane interposed between a lower and an upper compartment (CIM-plate for migration assay). In this way, the impedance-based detection of cell attachment, spreading and proliferation due to the gradual increase of electrode surface occupation may be monitored in real time and expressed as Cell Index. To determine the doubling time of A375, HT1080, and Sarc cell lines, cells re-suspended in growth medium were seeded on E-plates and impedance changes were continuously monitored for 70 h (Fig 1B). Only curves generated by seeding 4 × 103 cells/well were considered since those generated by seeding 2 × 103 cells/well did not reached a plateau until 90 h. According to their smaller sized dimension, A375 cells exhibited a long lasting adhesion/spreading phase and entered the growth phase (proliferation) and then stationary phase (plateau phase of growth) due to occupation of all entire microelectrode surface later, as compared to HT1080 and Sarc cells (Fig 1A and 1B). A375, HT1080 and Sarc cells reached the plateau phase after ≈70 h, 65 h, and 50 h, respectively (Fig 1B), and their doubling times calculated from the cell growth curve during the exponential growth were 32.8 ± 1.1 h, 16.2 ± 0.5 h and 10.87 ± 0.3 h, respectively (see S1 Fig). Since we did not employed cells subjected to cell cycle synchronization, we cannot exclude that, in the presence of serum, some cell division may occur on the bottom side of filter membranes, thereby affecting Cell Index. Therefore, to minimize the contribution of any cell division, cell migration experiments were performed for 12 h.

Fig 1

Experimental data.

A. Representative images of human melanoma A375, fibrosarcoma HT1080, or chondrosarcoma Sarc cells analysed by phase contrast microscopy. Original magnifications: 400x. Scale bar: 100 μm. B. Time-dependent proliferation of the considered human cell lines. Cells (2 × 103 cells/well) were seeded on E-plates and allowed to grow for 70 h in serum containing medium. The impedance value of each well was automatically monitored by the xCELLigence system and expressed as a Cell Index. Data represent mean ± SD (standard deviation) from a quadruplicate experiment. C. Cell migration of the indicated human cell lines monitored by the xCELLigence system. Cells were seeded on CIM-plates and allowed to migrate towards serum free medium (basal cell migration, black line) or medium plus 10% FBS. Cell migration was monitored in real-time for 12 h and expressed as Cell Index. Data represent mean ± SD from a quadruplicate experiment.

To evaluate cell motility in a system representative of the in vivo context, we compared the ability of A375, HT1080, and Sarc cell lines to migrate toward fetal bovine serum (FBS) which is a rich source of growth factor stimuli and chemotactic agents, which signal through binding to their cognate receptors [21]. To this end, cells (2 × 104 cells/well) were seeded on CIM-plates and allowed to migrate toward serum-free medium (basal cell migration) or growth medium, containing 10% FBS as a source of chemoattractants (directional migration), as described in [25]. As shown in Fig 1C, all cell lines exhibited a scarce basal cell motility (black lines), as their Cell Indexes did not change significantly along the time. On the other hand, all cell lines were able to respond to serum, although to a different extent. In agreement with their reported high motility [23, 24], both fibrosarcoma HT1080 and chondrosarcoma Sarc cells exhibited a comparable, high motility whereas a low response to FBS was retained by A375 cells (Fig 1C).

The mathematical model

In our mathematical model we schematize a single well of the CIM-plate used in the experiments as two cylindrical chambers, the upper and the lower chamber respectively, interfaced through the permeable membrane (Fig 2). The model considers two variables: the cell density and the amount of available FBS, that contains many chemotactic agents. Therefore, we have to take into account the possibility that some chemotactic agents may be degraded, consumed or internalized during the experiment. Thus, the FBS variable will describe the serum dynamics which includes a possible inactivation of some chemotactic agents. For both cells and the chemical signal we adopt a continuous description. This is justified by the high number of cells (or molecules) involved in the migration assay, typically in the order of 105 cell/cm3.

Fig 2

Schematic representation of a well of the CIM-plate.

An upper and a lower chamber are separated by a permeable membrane ΓM. In the migration assay in presence of chemoattractant, cells are placed in the upper chamber, and the chemoattractant is added in the lower chamber (directional migration). When measuring the basal migration experiment the well contains only cells (in the upper chamber) and a serum-free medium. In the mathematical formulation the spatial x−axis is oriented from the top to the bottom.

In the following, we first describe the rationale behind the proposed mathematical model, then we provide the explicit formulation in terms of equations. The cell population dynamics is the results of different contributions: diffusion, chemotaxis, spontaneous transport, cell adhesion/spreading. In particular, we consider a diffusion effect of cells in the environmental medium and the chemotactic effect of the FBS that attracts cells toward its higher concentrations. The transport term represents the so called basal migration, describing the cell transport through the pores of the permeable membrane, that we experimentally observe also in absence of chemotactic stimuli. The additional term of adhesion/spreading/proliferation represents the increase of Cell Index due to multiple reasons: better adhesion of cells to the biosensor, cell spreading, and possible cell proliferation. Such effects are indistinguishable, since the impedance-based estimate is related to the proportion of biosensor surface in contact with cells. Therefore, a better adherent, or spreaded cell, or duplicated cells produce an analogous increment of the surface contact. In this context, since we are limiting the observation time to 12 h (approximately or below the doubling time, see previous subsection), the observed effect can be related mostly to adhesion/spreading. The FBS variable is governed by a diffusion effect, coupled with a degradation term due to the chemotactic action (binding) during the cell migration. On the other hand an enzymatic degradation for FBS can be neglected in the considered experimental time range. The permeable membrane is modelled assigning the fluxes of the cell density and of the chemical signal through it, typically proportional to the difference of concentrations on the two sides of the interface.

Now we introduce the general equations of the mathematical model. Let Ω the domain consisting of the upper (ΩT) and lower (ΩB) chamber. We indicate with ΓT, ΓB and ΓM, respectively, the boundaries of the upper (top) chamber, of the lower (bottom) chamber, and the middle permeable membrane (see Fig 2). Then, let u(x, t) the cell density and φ(x, t) the FBS concentration, from the above considerations we have: where D, D, δ are positive constants, and χ(φ), g(u, φ) suitable functions, that will be specified later. On the boundary we assume the following conditions: where n is the downward normal versor, k is a constant, while k(u) depends on the cell density, and finally uT, φT, uB, φB are the limit values of u and φ on the interface ΓM from the upper and lower chamber, respectively.

Let us now specify the contribution of the different terms in the first and second equation of the proposed system, Eq (1)1 and Eq (1)2. For the diffusion term in Eq (1)1 we assume a constant diffusion coefficient D. The chemotaxis term involves a modulating function χ(φ), which takes into account a possible saturation effect for high concentration of chemoattractant. A possible choice is with χ1 and χ2 positive constants. Similar modelling functions can be found, for example, in [26]. The spontaneous transport of the cell in absence of chemoattractant is modelled as a transport term at velocity Vtransp. We can assume a constant velocity in the direction of the vector n as the limit velocity achieved by the cells in the viscous environment. About the cell adhesion/spreading effect, we consider the function that establishes a logistic growth for u, as it can be deduced from related experiments of proliferation (see previous subsection, Fig 1, and section Materials and methods). The function W(x) is a weight function, which spatially forces the spreading effect as described in the following. Firstly, from the experimental point of view we observe that when cells migrate in the lower chamber, they remain adherent to the bottom side of the membrane. However, for simplicity reasons, in the proposed model cells crossing the membrane are not confined on its lower surface, but they can freely move in the lower chamber. Therefore, to obtain the number of migrated cells, it is necessary to consider the entire lower well, integrating the cell density on it. In our framework, the adhesion/spreading effect involves the cells on the upper face of the membrane and also all cells crossed into the lower chamber. Therefore, a possible choice for the W function, along the x-axis, is where xM is the position of the central membrane and a suitable constant. In the following we will assume in the order of two cell diameters. The term in Eq (7) considers that FBS serum promotes the cell adhesion/spreading through its growth factors (previous subsection), possibly with a saturation effect, while in its absence (for example in the basal migration experiment) such increase in cell density is assumed negligible. The constants α1, α3 can be estimated, for a specific cell line, fitting related proliferation data obtained at concentration of FBS . The underlining assumption of these estimates is that proliferation assays and migration assay show similar adhesion/spreading rate at least in the earlier times (see next subsection, and section Materials and methods). We recall that if the time of observation of the migration assay remains limited in the interval of 12 h, the increment of Cell Index, can be mostly attributed to the adhesion/spreading effect. However, on higher times the contribution of Eq (7) could be able to reproduce also an increase in cell density due to cell proliferation (see also section Discussion and conclusions).

For the FBS signal in Eq (1)2 we have a diffusion term with constant coefficient D. Along with this, we consider a serum consuming term proportional to the product between cell density and chemical signal, which takes into account the inactivation of the chemotactic agents due to binding process. Conversely, on the scale of the examined experiment the molecular degradation of the serum can be neglected.

About the boundary conditions, Eqs (2) and (3) represent zero flux for the cell and for the chemoattractant on the top and bottom side of the well, since no mass leaves our domain. In Eqs (4) and (5) we fix Kedem-Katchalsky boundary conditions, meaning that the flux of cells and FBS through the membrane is proportional to the difference between the concentrations at the top (uT) and bottom (uB) sides of the boundary (see [27], and [28] for a mathematical and modellistic treatment of this conditions). For the φ signal we assume a constant transmission coefficient, while for u the transmission term is considered as a function of the cell density. In particular in k(u) we assume possible crowding effects on both sides of the interface. A suitable function can be where k, k, k are constants. Notice that Eq (9) decreases for increasing cell density on the membrane. In particular, we have two contributions in the denominator: one given by the cell density on the upper side of the interface Γ (uT); the other given by a similar contribution on the lower side of Γ, possibly up to the power p. As we have observed above, we need to integrate the cell density on the entire lower chamber ΩB. Numerical data suggest p = 2 as a suitable power, which we will assume in the following. All the above considerations are then summarized in the following system of equations: Initial concentrations for u(x, t) and φ(x, t) will be in the form Ω being the portion of the upper chamber, with positive cell density at t = 0.

Parameter estimation and sensitivity analysis of the mathematical model

In our numerical tests we applied the mathematical model to three different cell lines: Sarc, HT1080, A375; and two conditions: migration toward chemoattractant (FBS in our case) and basal migration. The second condition corresponds to choose χ1 = α1 = 0 in the system (10), so that Eq (10)1 and Eq (10)2 decouple, and we can simulate only the dynamics of the cell density.

For symmetry reasons, we can simulate a one-dimensional version of system (10), lengthwise the cylindrical domain. Such assumption is in agreement to the impedance-based measurement of the Cell Index performed by the cell analyser, and used to compare numerical data.

In order to simulate the dynamic of the model, all parameters have to be chosen. To this purpose we remark that some of them are already available in biological or modellistic literature, while the others have been calibrated on the experimental data. Table 1 summarizes the set of parameters used in our simulations for the different cell lines. For those retrieved from scientific papers we provide the reference in the last column, while for the others, marked as “data driven”, we specify the experiment (i.e. proliferation, migration, or basal migration) from which we have derived them.

Table 1

Initial data and parameters of the mathematical model.

Initial datum or parameter	Definition	Estimated value	Source
u0	initial maximum cell denity in Eq (11)	≈30200, ≈45300, ≈60400 cellcm−1 (Sarc) ≈30200 cellcm−1 (HT1080, A375)	Exp. setup: sec. Mat. and Meth.
φ0	initial maximum FBS concentration in Eq (12)	18.39 μlcm−1	Exp. setup: sec. Mat. and Meth.
Du	cell diffusion	1 × 10−3 cm2 h−1 (Sarc) 2.5 × 10−3 cm2 h−1 (HT1080) 8 × 10−4 cm2 h−1 (A375)	data driven from basal migr. exp.
Dφ	FBS diffusion	3.7 × 10−3 cm2 h−1	[29]
χ1	first chemotactic constant	3 × 10−3 cm3 μl−1 h−1 (Sarc) 2.5 × 10−3 cm3 μl−1 h−1 (HT1080) 1 × 10−3 cm3 μl−1 h−1 (A375)	data driven from migr. exp.
χ2	second chemotactic constant	4.75 × 10−8 μlcm−1 (Sarc, HT1080, A375)	data driven from migr. exp.
Vtransp	transport velocity	9 × 10−3 cmh−1 (Sarc) 2 × 10−3 cmh−1 (HT1080) 1.3 × 10−9 cmh−1 (A375)	data driven from basal migr. exp.
α1	logistic growth coefficient	0.154 h−1 (Sarc) 0.135 h−1 (HT1080) 0.118 h−1 (A375)	data driven from prolif. exp.
α2	dependence on FBS in the logistic growth	10−6 μlcm−1 (Sarc, HT1080, A375)	data driven from migr. exp.
α3	limit value in the logistic growth	2.08 × 105 cellcm−1 (Sarc) 2.26 × 105 cellcm−1 (HT1080) 1.04 × 105 cellcm−1 (A375)	data driven from prolif. exp.
φ¯	FBS concentration in proliferation experiments	19.64 μlcm−1	Exp. setup: sec. Mat. and Meth.
δ	FBS degradation	10−8 cmh−1 cell−1 (Sarc, HT1080) 3.5 × 10−5 cmh−1 cell−1 (A375)	data driven from migr. exp.
ku1	cell transmission coefficient on the membrane	2 cmh−1 (Sarc, HT1080, A375)	data driven from basal migr. exp.
ku2	crowding coefficient on the upper side of the membrane	1 × 10−5 cmcell−1 (Sarc) 5 × 10−5 cmcell−1 (HT1080) 2 × 10−8 cmcell−1 (A375)	data driven from basal migr. exp.
ku3	crowding coefficient on the lower side of the membrane	6 × 10−8 cell−2 (Sarc, HT1080, A375)	data driven from basal migr. exp.
kφ	FBS transmission coefficient on the membrane	8.8 × 10−2 cmh−1 (Sarc, HT1080, A375)	data driven from migr. exp.

Estimates of initial data, physical and biological parameters. About the model parameters, values were retrieved from scientific literature, or estimated from proliferation or migration assays. For those obtained from migration experiments we used 2 × 104 cells/well.

In detail, the constants u0, φ0, were assigned by the experimental protocol (section Materials and methods). The coefficient D was set according to [29]. Coefficients related to the cell proliferation, i.e. α1, α3, were obtained, for a specific cell line, from proliferation experiments. This kind of assays was performed in real time on E-plates, using the same technology of the migration ones, for each cell line and at a known FBS concentration (section Materials and methods). Experimental curves showed a logistic growth in the cell density, and were interpolated to estimate the above mentioned parameters of our interest. Then, the parameters which do not involve chemotactic or growth effects, that are D, Vtransp, k, k, k, were calibrated on the basal migration curves, fixing in the model χ1 = α1 = 0. Finally, χ1, χ2, α2, δ, k, were calibrated consistently with the other parameters, on the migration curves in presence of chemoattractant.

We observe that, as in many mathematical models of biological phenomena, the lack of complete information from the experiments on the parameter values necessarily imposes an uncertainty in the response of the model. To obtain as reliable results as possible, we have studied the influence of the parameters on the behaviour of the model through a local sensitivity analysis [30, 31], as described below. Such approach allows us to estimate an influence index between the variation of a parameter and a particular observed output of the model. In our analysis we consider the variation of a single parameter at a time, so interactions among coefficients are neglected. This is useful for a first exploration of the parameter space.

Let p0 a parameter value and ε a small deviation on p0, let f(p0) an output obtained for the p0 value, we defined the sensitivity index S as the following ratio between relative variations:

Table 2 shows the S value in Eq (13) for the parameters that we calibrated on the experimental data. The small deviation ε was assumed equal to 0.05p0, that is a 5% deviation on the parameter value. The observed output f was the Cell Index at the final time of observation, corresponding to 12 h. Moreover, Table 2 shows also, in the second column, the percentage variation of the examined parameter given by

Table 2

Sensitivity analysis for the parameters of the mathematical model.

Parameter variation	Cell Index variation at 12 h	S
Du + ε	+0.38%	0.08
Du − ε	−0.39%	0.08
χ1 + ε	+0.90%	0.18
χ1 − εχ2 + ε	−0.97%−1.2 × 10−4%	0.192.3 × 10−5
χ2 − ε	+1.2 × 10−4%	2.3 × 10−5
Vtransp + ε	+0.59%	0.12
Vtransp − ε	−0.59%	0.12
α2 + ε	< 10−5%	< 10−5
α2 − ε	< 10−5%	< 10−5
δ + ε	+7.9 × 10−3%	1.6 × 10−3
δ − ε	−8 × 10−3%	1.6 × 10 − 3
ku1 + ε	+5.2 × 10−2%	1 × 10−3
ku1 − ε	−5.7 × 10−2%	1.1 × 10−3
ku2 + ε	−7.5 × 10−5%	1.5 × 10−5
ku2 − ε	+7.5 × 10−5%	1.5 × 10−5
ku3 + ε	−5.2 × 10−2%	1 × 10−2
ku3 − ε	5.3 × 10−2%	1 × 10−2
kφ + ε	+0.15%	0.03
kφ−ε	−0.17%	0.03

Local sensitivity analysis for parameters in Table 1 calibrated from numerical simulations. Second column shows the relative percentage variation as in Eq (14), choosing as observed output f the Cell Index at the final time of the simulation (12 h), and considering ε corresponding to a 5% variation. Third column contains S in Eq (13).

Numerical simulations on basal and directional cell migration of three different cell lines

In this section we show the performance of our dynamical model in describing experimental results. In order to compare our numerical data with those obtained by xCELLigence analyser we express the cell density in term of Cell Index. The Cell Index linearly depends on the cell density [32], and the coefficients of this linear dependence can be estimated from an experiment of cell proliferation, in which a known number of cells is placed on an impedance-based biosensor and their Cell Index measured in time (section Materials and methods).

In our simulations we chose the domain Ω = [0, 1.8] (cm) according to the well height in the used CIM-plate, which permeable membrane is placed in the middle, at x = 0.9 cm [33]. As observed in previous sections, the observation time was fixed at 12 h (Fig 1 in section Results). For the space discretisation, to preserve stability we adopted a non-uniform mesh. In particular, we fixed Δx = 10−2 cm, while in proximity of the membrane we reduced the spatial step to the finer Δxf = 10−6 cm. The time step was chosen as the maximum value able to ensure stability and non-negativity of the solution, that is Δt = 10−3 h (for further details see section Materials and methods).

In all numerical tests, the parameters of the model, estimated as described in previous section, were chosen according to Table 1. Dynamical simulations were compared with the relative experimental curves computed as described below. For each cell line at least three independent experiments were available. Each experiment was performed in quadruplicate on the same CIM-plate, and the xCELLigence data were recorded as mean value (section Materials and methods). We consider as resulting experimental curve for each cell line, the average of the independent replicates (see S2 Fig for full raw data).

For each comparison we estimated also the relative MSE error, given by where n is the number of experimental time steps, and , c are respectively the numerical and the experimental Cell Index. When necessary, was interpolated on time steps of c. In the following we will indicate with MSEmigr and MSEbasal respectively the MSE relative to the migration and basal migration simulations.

Fig 3 shows a numerical simulation of the model (10) with parameters fixed as in Table 1 in comparison with the experimental data. Specifically, panel (a) and (b) refer to Sarc cell line, reporting results for basal migration (a) and full system (10) (b) respectively. Experimental curves are marked in red for basal migration, and green for chemotactic migration. Both panels display the Cell Index curve versus time. The estimate of the relative MSEs are given by MSEbasal = 0.0376 and MSEmigr = 0.0052. Fig 3(c) and 3(d) refer to HT1080 cell line. In this case we obtain the values MSEbasal = 0.0166 and MSEmigr = 0.0068. Finally, in Fig 3(e) and 3(f) we consider the A375 cell line. For the relative MSE we estimate MSEbasal = 0.0083 and MSEmigr = 0.0054.

Fig 3

Numerical simulations on Sarc, HT1080, and A375 cell lines.

For each cell line, panels on the left (a),(c),(e) show the basal migration in absence of chemoattractant. Numerical curves (blue) were compared with experimental data (red). Panels on the right (b),(d),(f) show the migration curves. The simulated values of Cell Index (blue) were compared with experiments (green). Here and in the following figures the experimental curves were obtained as the average of at least three experiments in quadruplicate (S2 Fig). About the MSE value on the Cell Index, defined in Eq (15), we estimated, respectively, the following values: panels (a)-(b) MSEbasal = 0.0376 and MSEmigr = 0.0052; panels (c)-(d) MSEbasal = 0.0166 and MSEmigr = 0.0068; panels (e)-(f) MSEbasal = 0.0083 and MSEmigr = 0.0054.

Confirming the mathematical model on chondrosarcoma Sarc cells

In previous section we have shown that, after a suitable parameter calibration, the proposed mathematical model was able to describe the cell migration of three different cell lines, with a very good concordance with the experimental data. Here we investigated the model capability to make predictions about new experiments. To this aim, we used our model to predict the Cell Index on Sarc cell lines performed with different numerosities of cells. Therefore, we applied our mathematical model, using the parameters estimated on the Sarc cell line in the case of 2 × 104 cells in migration (Table 1), and we estimate the behaviour corresponding to 3 × 104 and 4 × 104 cells/well. Then, in related experiments, cells were seeded at these two different densities on CIM-plates and allowed to migrate towards serum-free medium (basal cell migration) or medium plus 10% FBS. Cell migration was monitored in real-time for 12 h as changes in Cell Index. In Fig 4 we show the comparison between these numerical curves and Cell Index data obtained by the xCELLigence analyser in the case of migration towards chemoattractant. The displayed experimental data represent an average value of three and four different experiments, respectively for the case of 3 × 104 and 4 × 104 cells/well (see S2 Fig). In all cases we found a nice agreement with the experimental evidences. In particular, for 3 × 104 and 4 × 104 cells/well, we estimated respectively the relative MSE value as MSEmigr = 0.0077, and MSEmigr = 0.0183.

Fig 4

Sarc chondrosarcoma cell line. Confirming the mathematical model.

Model (10) was simulated with parameters fixed as in Table 1, obtained with 2 × 104 cells/well, and varying the initial cell density u0. In (a) and (b) numerical data of migration curves were compared with experimental Cell Index respectively in the case of 3 × 104 and 4 × 104 initial cell number. MSE value on the Cell Index was estimated in MSEmigr = 0.0077 and MSEmigr = 0.0183, respectively in (a) and (b).

Discussion and conclusions

Cell migration is a process that offers rich targets for intervention in key pathologic conditions, including cancer. Indeed, the development of metastases requires the activation of a series of physiological and biochemical processes that govern the migration of tumour cells from the primary tumour site, the invasion through the basement membrane, the entry of metastatic cells into the blood vessels and finally localization to the second site [1]. Therefore, targeting cell motility has been increasingly accepted as a new approach for the clinical management of metastatic patients and in the future, quantitative analysis of the motility of tumour cells derived from cancer patients could provide a new potential parameter predictive of patient outcomes. The recent expansion of mathematical modelling is already contributing to cancer research by helping to elucidate mechanisms of tumour initiation, progression and metastases as well as intra-tumour heterogeneity, treatment responses and resistance [12]. Parametrization of cell motility is often difficult given the available experimental model systems. With the advent of high throughput systems, there has been a movement towards the use of a number of cell-based assays useful for studying cell migration. A recent technology, xCELLigence RTCA, has been increasingly accepted as a platform for high-throughput determination of cell motility dynamics in real time using micro-electronic biosensors [34]. In this paper we propose a macroscopic mathematical model, based on convection-reaction-diffusion equations, for the cell migration assay.

Previous mathematical models on in vitro Boyden-like assays dealt mainly with the invasion experiment [19-21]. Among these, the authors in [21] studied cancer cell invasion through a theoretical model compared with real-time impedance-based assays. By contrast, in this paper, we proposed a PDEs model in relation to the cell migration experiment relying on the xCELLigence real-time technology. Our model differs from [21], where it was assumed that the simulated Cell Index is proportional to the fraction of cells that reaches the upper well bottom. The authors assumes that the pores dimension of the permeable membrane are larger enough to allow cells, quite easily, to cross it. On the contrary, cell lines employed in our migration experiments present much larger dimensions than membrane pores (8 μm) (Fig 1A). For this we consider the effective cell crossing through the porous interface, and the simulated Cell Index was computed on the basis of the fraction of cells migrating in the lower chamber. Moreover, our transmission coefficient in the boundary conditions includes also possible crowding effects, being assumed as a decreasing function of the cell density on both the faces of the separating membrane. Finally, to model the basal migration effect, as described in section Results, our model considered also a spontaneous transport of cells across the permeable membrane, present even in absence of chemotactic stimuli, that is not considered in [21]. This allowed us to recover experimental data through the basal experiment and to estimate on it some parameters to be included in the full migration model.

Numerical simulations has been performed to compare the model dynamics with experimental raw data obtained by the xCELLigence RTCA in absence or presence of a chemotactic gradient. We also validate the performance of our model by comparing the results of simulations with other experimental data, on chondrosarcoma Sarc cell line, not used for estimating model parameters. Numerical findings showed a nice agreement with the acquired experimental data. Therefore, overall we can infer that tumour cells migration can be described using mathematical models as a predictable process dependent on biophysical laws and experimental parameters.

Starting from the present paper, some interesting issues can be investigated as future perspectives. Firstly, we could explore the quantitative and qualitative accuracy of the model to simulate different experimental conditions in a migration assay, such as the initial serum concentration, or to test the effects of various chemoattractants. In this regard, it could also be interesting to introduce the action of chemotactic inhibitors on the cell motility. In the future, we will explore the possibility of simulating in silico the ability of inhibitors of cell migration to counteract the motility of primary tumour cells derived from patients affected by solid tumours, in order to design more personalized therapeutic strategies.

Materials and methods

Cell Lines

Human melanoma A375 cell line purchased from American Type Culture Collection (ATCC) was cultured in RPMI 1640 medium (Lonza, Milan, Italy), supplemented with 3 mM L-glutamine (Invitrogen-Gibco®/Life Technologies, Monza, Italy), 2% penicillin/streptomycin and 10% fetal bovine serum (FBS). Highly mobile human fibrosarcoma HT1080 cell line [23], also purchased from ATCC, and human chondrosarcoma Sarc cells derived from a chondrosarcoma primary culture [24], preliminary characterized for their ability to migrate toward several chemoattractants [35], were cultured in Dulbecco Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS), 100 IU/ml penicillin and 50 μg/ml streptomycin. All cells were maintained at 37°C in a humidified atmosphere of 5% CO2.

Cell Proliferation

Cell proliferation was assessed using the xCELLigence RTCA technology as described [36]. For these experiments, the impedance-based detection of cell attachment, spreading and proliferation was assessed by using E-plates which are provided of microelectrodes attached at the bottom of each well. First, 100 μl of growth medium was added to each well, the plate was locked at 37°C in a humidified atmosphere of 5% CO2 and the background impedance was measured. Cells were counted, suspended in in 100 μl growth medium, seeded (2 × 103 or 4 × 103 cells/well) and allowed to grow for 70 h. The impedance value of each well was automatically monitored by the xCELLigence system and expressed as a Cell Index.

Cell Migration

Cell migration was monitored using the xCELLigence RTCA technology as described in [36]. For these experiments, the impedance-based detection of cell migration was assessed using CIM-plates which are provided of interdigitated gold microelectrodes on bottom side of a microporous membrane (containing randomly distributed 8 μm-pores) interposed between a lower and an upper compartment. Briefly, 160 μl of serum-free medium with/without 10% FBS and 30 μl of serum-free medium were added to the lower and upper chambers, respectively, prior to lock the plate at 37°C in a humidified atmosphere of 5% CO2 for 60 minutes (to obtain the equilibrium between the two compartments), according to the manufacturer’s guidelines. Then, background signals generated by cell-free media were measured, detached cells were counted, suspended in 100 μl serum-free medium and seeded (2 × 104, 3 × 104, 4 × 104 cells/well) in the upper chamber. Microelectrodes detect impedance changes which are proportional to the number of migrating cells and are expressed as Cell Index. Cell migration was monitored in real-time for 12 h. Each experiment was performed at least three times in quadruplicate. Raw data are given in S1 and S2 Files.

Numerical methods

The numerical approximation scheme used in the simulation of the model (10) employed a finite difference method on a spatial domain Ω = [a, b], consisting of upper and lower domains ΩT, ΩB, interfaced through the membrane ΓM. Let Δx, Δt the space and time steps, we defined the grid points (x, t), where x = iΔx and t = kΔt. The approximation of a function f(x, t) at the grid point (x, t) was denoted as . To ensure non-negativity in the numerical simulations, due to the boundary conditions on the permeable membrane, in its proximity we needed to discretise our equations on a finer spatial mesh x = jΔxf, Δxf < Δx.

For the diffusion Eq (10)2 we applied on the internal nodes a central scheme in space and an implicit scheme in time for the diffusive term, while the reaction term was put in explicit: Similarly on the finer mesh with spatial step Δxf.

For the advection-diffusion Eq (10)1, let we assumed and for the internal nodes we adopted the scheme with function W(x) defined in Eq (8), and where the last term introduced an artificial viscosity in order to preserve scheme stability (see for example [37]).

For the boundary conditions (10)3,4 on ΓT (x = x0) and on ΓB (x = x) we used the second order one-sided approximation of the second derivative in the form: On the boundary ΓM (x = xM), let uT, uB the variable u in ΩT, ΩB respectively. From Eq (10)5,6, for u we employed where was given by Eq (18) and Similarly, Eq (10)6 was discretised as In our numerical simulations we used Δx = 10−2 cm, while the interval [xM − Δx, xM + Δx] centred on the membrane was discretised with Δxf = 10−6 cm. Stability and non-negativity of numerical solutions were obtained by choosing Δt = 10−3 h.

Doubling times of Sarc, HT1080, and A375 cell lines.

Cells (2 × 103 cells/well) were seeded on E-plates and allowed to grow for 70 h in serum containing medium. The impedance value of each well was automatically monitored by the xCELLigence system and expressed as a Cell Index. Doubling times were calculated, using the xCELLigence RTCA software, from the cell growth curves during exponential growth given in round brackets for each cell line. Doubling time is expressed in term of mean value ± SD (standard deviation) from a quadruplicate experiment.

(TIF)

Click here for additional data file.

Cell Index data recorded by xCELLigence of the different experiments in our study.

Panels (a),(c),(e),(g),(i) describe the basal migration, (b),(d),(f),(h),(j) the migration in presence of FBS. In each panel the curves represent an independent experiment carried out in quadruplicated and averaged. The observed curves in Figs 3 and 4 are obtained as the average of the curves showed here.

(TIF)

Click here for additional data file.

Basal migration xCELLigence raw data.

The file contains 19 different spreadsheets organized with respect to cell lines, initial cell numbers, and independent experimental replicates. Within the same spreadsheet the first column contains the time (in hours), second and third column contain the mean basal migration Cell Index of a quadruplicate experiment and its standard deviation, respectively.

(XLSX)

Click here for additional data file.

Migration xCELLigence raw data.

The file contains 17 different spreadsheets organized with respect to cell lines, initial cell numbers, and independent experimental replicates. Within the same spreadsheet the first column contains the time (in hours), second and third column contain the mean migration Cell Index of a quadruplicate experiment and its standard deviation, respectively.

(XLSX)

Click here for additional data file.