Multiscale modeling of collective cell migration elucidates the mechanism underlying tumor-stromal interactions in different spatiotemporal scales.

Metastasis is the pathogenic spread of cancer cells from a primary tumor to a secondary site which happens at the late stages of cancer. It is caused by a variety of biological, chemical, and physical processes, such as molecular interactions, intercellular communications, and tissue-level activities. Complex interactions of cancer cells with their microenvironment components such as cancer associated fibroblasts (CAFs) and extracellular matrix (ECM) cause them to adopt an invasive phenotype that promotes tumor growth and migration. This paper presents a multiscale model for integrating a wide range of time and space interactions at the molecular, cellular, and tissue levels in a three-dimensional domain. The modeling procedure starts with presenting nonlinear dynamics of cancer cells and CAFs using ordinary differential equations based on TGFβ, CXCL12, and LIF signaling pathways. Unknown kinetic parameters in these models are estimated using hybrid unscented Kalman filter and the models are validated using experimental data. Then, the principal role of CAFs on metastasis is revealed by spatial-temporal modeling of circulating signals throughout the TME. At this stage, the model has evolved into a coupled ODE-PDE system that is capable of determining cancer cells' status in one of the quiescent, proliferating or migratory conditions due to certain metastasis factors and ECM characteristics. At the tissue level, we consider a force-based framework to model the cancer cell proliferation and migration as the final step towards cancer cell metastasis. The ability of the multiscale model to depict cancer cells' behavior in different levels of modeling is confirmed by comparing its outputs with the results of RT PCR and wound scratch assay techniques. Performance evaluation of the model indicates that the proposed multiscale model can pave the way for improving the efficiency of therapeutic methods in metastasis prevention.

Introduction

Cancer cell metastasis is a key step of cancer progression that accounts for the majority of patient deaths[1,2]. Despite the inherent complexity of metastasis stages, recent studies show that not only intracellular reactions in cancer cells, but also mutual interactions between cancer cells and tumor microenvironment (TME) orchestrate the events during a metastasis cascade[3,4]. In other words, studying the behavior and characteristics of cancer cells without considering their interactions with TME components such as extracellular matrix (ECM), blood vessels, immune cells, adipose cells, and fibroblasts does not provide comprehensive information about the causes and circumstances of tumor cells movement[5,6].

Cancer associated fibroblasts (CAFs) as one of the abundant factors in TME have a major impact on tumor behavior and they are known as principal participants in tumor growth and invasion[7,8]. They directly contribute to cancer cell proliferation, tumor growth, and invasion through stimulation of various growth factors and chemokines such as transforming growth factor β (TGFβ), vascular endothelial growth factor (VEGF), leukemia inhibitory factor (LIF), and C-X-C motif chemokine 12 (CXCL12) in context-dependent manner[8-10]. These factors are responsible for initiation of signaling pathways and consequent biochemical reactions in both cancer cells and CAFs -as a principle provocative components of tumor metastasis- which lead to tumor progression and metastasis. To address the effects of CAFs on tumor metastasis, it is necessary to investigate the communications between CAFs and cancer cells, as well as their intercellular interactions which result in tumor invasion.

Cancer metastasis is a complex multiscale process in which reactions occur in different time and space scales[11]. Furthermore, the number of these interactions are relatively high and the related signaling pathways usually have crosstalk with each other. Therefore, an advanced tool is required not only for detailed study of these interactions, but also for efficient analysis of the involved mechanisms in different stages of metastasis. This is where systems biology comes into play and provides useful techniques for investigating the complex biological behaviors such as mathematical modeling approaches. Using these methods, we obtain the opportunity to understand the behavioral kinetics inside the cells as well as communications between different cell types in TME. Moreover, mathematical modeling tools enable us to explain how interactions through three different scales including molecular, cellular, and population lead to tumor migration.

Over the last decades, various mathematical modeling approaches have been applied for revealing tumor initiation and growth mechanisms[12-17]. However, the number of studies that model metastasis based on communication of tumor and its milieu, are relatively low. In[18], a mathematical model has been presented for describing cancer cell intravasation. Using a multiscale approach and a hypothetical framework for biological shape of blood vessel, the model represents the effect of cadherin protein pathways in cancer cell migration.

Also, the role of ECM remodeling in cancer cell invasion by means of degradation enzymes has been investigated in some recent studies[19-21]. Kim et al.[22] developed a hybrid model that incorporates the interactions of stromal and cancer cells in TME to show how theses communications play role in early development of tumor and its invasion to the stroma. In addition, a number of mathematical models based on agent-based methods have been applied for determining the interactions of platelet cells, stem cells, and tumor cells, as well as the role of immune cells in cancer metastasis[23-26].

It should be noted that there are a limited number of studies based on mathematical modeling that focus on the effects of CAFs on cancer cells’ behavior. For example, based on the singular value decomposition approach, regulatory rules for interactions of cancer cells and adjacent fibroblasts were modeled in[27]. In addition, a mathematical model for the transformation of fibroblasts into myofibroblasts in TME, as well as the effect of myofibroblasts on tumor cell proliferation via Epidermal growth factor (EGF) release, has been presented[28]. Furthermore, in a recent study we presented an agent-based model of fibroblast switching behavior in tumor suppression and progression in normal and active states[29]. Nonetheless, the mutual conversation between CAFs and tumor cells at microscopic and macroscopic levels that leads to metastasis has never been mathematically investigated.

In the current study, we present a multiscale model of cancer cell metastasis due to their interactions with TME components such as CAFs and ECM. At first, a mathematical modeling framework has been developed to represent the dynamics of cancer cells and CAFs individually based on ordinary differential equations (ODEs). The dynamical models represent biochemical reactions in TGFβ, LIF, and CXCL12 signaling pathways which were retrieved from the literature. As a result, each cell type has its own set of ODEs at the molecular level. Similar to our previous work in[29], hybrid unscented Kalman filter (HUKF) is used to unravel the unknown parameters as a reliable parameter estimation approach[30]. Then, using partial differential equations (PDEs) for describing the spatial distribution of molecules throughout the TME, we extended the ODE systems to a coupled ODE–PDE system. This coupled system can explain the communications between cancer cells, CAFs and ECM in a three-dimensional space. The proposed model can explain the nonlinear effect of CAFs on the invasive activity of cancer cells at the cellular level.

Moreover, the coupled ODE–PDE system can designate the state of each individual cancer cell as quiescent, proliferative, or migratory. The states of each cell are determined due to the concentration level of some significant factors in signaling pathways such as CXCL12, and SNAIL. Finally, we developed a multiscale model which illustrates the movement of cancer cells in a three dimensions domain using a coupled ODE–PDE system and a biomechanical model at population level. This model represents the movement of cells due to the forces applied to them. Finally, the model is validated in both molecular and population levels by comparing the simulation results with the data extracted from RT PCR and migration assay. The proposed nonlinear multiscale model can not only be utilized in determining the role of CAFs in metastasis, but also might pave the way to finding advanced cancer treatment methods.

Material and methods

Multiscale model of cancer metastasis

The proposed model in this study is comprised of multiple layers each representing interactions at the molecular (within a cell), cellular (between cells), and population levels. The multi-layer model of cancer cell metastasis includes interactions in different components such as cancer cells, CAFs, and their mutual interactions in the TME. At the molecular level, nonlinear ODE models are used to represent the dynamics of TME components which are based on reactions in some signal transduction pathways. Then, the cellular level is basically characterized by the signaling molecules between cancer cells and CAFs. These signals show the impacts of cancer cells on CAF functioning and consequently, the effects of CAFs on cancer cell invasion. In this level, the dynamics show how the behavior of TME components are controlled under the influence of each other. Finally, cancer cell migration is depicted in three dimensions space using a Newtonian dynamic at the population level of the model. In order to model the movement of cancer cells, different layers are connected to each other appropriately. The various levels of modeling technique are given below. Figure 1 depicts the various levels of our multiscale model as well as the actions required at each level to model cancer cell metastasis.

Figure 1

Flowchart of necessary steps in each level of multiscale model of cancer cell metastasis.

Dynamical modeling of intracellular interactions within cancer cells and CAFs

At the molecular level, we present two unique nonlinear ODE models for characterizing the dynamics of cancer cell and CAF as two different cell types in TME. Unlike normal cells, some signaling pathways which are responsible for cell growth, proliferation, and motion function abnormally in tumor cells[31]. As a result, we can obtain abundant information from the well-known growth factor expression variations and their linked signaling pathways for better understanding cancer mechanisms. Due to vital roles of TGFβ, LIF and CXCL12 signal transduction pathways in tumor growth, proliferation, motility, and invasion, the ODE models are based on biochemical reactions involved in these signaling pathways[32-35]. They are also the key activators of fibroblasts in TME and abundant within both cancer cells and CAFs[36,37].

As demonstrated in Figs. 2 and 3, dynamical models are mathematical representations of biochemical reactions in TGFβ and CXCL12 pathways, as well as TGFβ and LIF signaling cascades with their crosstalk for describing the dynamics of cancer cells and CAFs, respectively. TGFβ, as one of the important factors released in the TME by both CAFs and cancer cells, binds to its receptors (TGFβRs) and consequently initiates a chain of reactions inside the cell. The constitution of active ligand-receptor compounds leads to phosphorylation of SMAD2/3 factors. The phosphorylated SMAD2/3 and SMAD4 then link together to create a heterodimeric complex. This complex imports to the nucleus and acts as a transcription factor for downstream genes such as SMAD7 and Zinc finger protein SNAI1 (SNAIL). SMAD7 inhibits SMAD2/3 phosphorylation upstream in the pathway, resulting in a negative regulatory feedback in TGFβ signaling[38,39]. It should be noted that LIF, C-X-C chemokine receptor type 4 (CXCR4), CXCL12, and TGFβ itself are expressed as the target genes of the pathway[32,33,40,41]. Using this fact, we can conclude that TGFβ interferes in LIF and CXCL12 signaling pathways by increasing the level of these factors, and also it can drive its pathway in both paracrine and autocrine fashions[32,42]. These two points are reflected in our modeling procedure.

Figure 2

Biochemical reactions involved in dynamical model of TGFβ and CXCL12 pathways as well as their crosstalk within a cancer cell. The name of each reaction is shown on the corresponding arrow. More information about the mathematical representation of reactions are included in Tables 1, S1, and S3.

Figure 3

Biochemical reactions involved in dynamic model of TGFβ and LIF pathways as well as their crosstalk within the CAFs. The name of each reaction has been shown on the corresponding arrow. More information about the mathematical representation of reactions are included in Tables 2, S2, and S4.

Table 1

Biochemical reactions and their mathematical representations in cancer cell dynamics.

Description of the interaction	Mathematical expression
Production of TGFβ receptor	v_1=k1+x_1
Degradation of TGFβ receptor	v_2=k1-x_1
Association of TGFβ-TGFβR complex	v_3=k2+xu_3x12
Dissociation of TGFβ-receptor complex	v_4=k2-x_2
Production of cytoplasmic SMAD3	v_5=k3+x_3
Degradation of cytoplasmic SMAD3	v_6=k3-x_3
Production of cytoplasmic SMAD4	v_7=k4+x_5
Degradation of cytoplasmic SMAD4	v_8=k4-x_5
Phosphorylation of SMAD3	v_9=k5+x_2x_3x_3+K_s1
Dephosphorylation of SMAD3	v_10=k5-x_4
Association of pSMAD3-4 Complex	v_11=k6+x_4x_5
Dissociation of pSMAD3-4 Complex	v_12=k6-x_6
Nuclear import of pSMAD3-4 Complex	v_13=k_i7x_6
Nuclear export of pSMAD3-4 Complex	v_14=k_e7x_7
Degradation of pSMAD3-4 Complex	v_15=k7-x_7
Production of SMAD7 in the pathway	v_16=k8+x_7
Degradation of SMAD7	v_17=k8-x_8
Inhibitory effect of SMAD7 on pSMAD3	v_18=k8i-x_8
Production of TGFβ in the pathway	v_19=k9+x_7
Degradation of TGFβ	v_20=k9-xu_1
Association of TGFβ-TGFβR complex	v3^′=k2+xu_1x12
Production of LIF in the pathway	v_21=k10+x_7
Degradation of LIF	v_22=k10-xu_2
Production of CXCR4 in the pathway	v_23=k11+x_7
Degradation of CXCR4	v_24=k11-x_9

Table 2

Biochemical reactions and their mathematical representations in CAF dynamics.

Description of the interaction	Mathematical expression
Production of CXCL12 in the pathway	v_23=k11+x_7
Degradation of CXCL12	v_24=k11-xu_4
Production of LIF receptor	w_1=h1+x_10
Degradation of LIF receptor	w_2=h1-x_10
Association of LIF-LIFreceptor complex	w_3=h2+x_10xu_2
Association of LIFcaf-LIFreceptor complex	w_3=h2+x_10x_9
Dissociation of LIF-LIFreceptor complex	w_4=h2-x_11
Production of JAK	w_5=h3+x_12
Degradation of JAK	w_6=h3-x_12
Phosphorylation of JAK by LIF-receptor complex	w_7=h4+x_11x_12x_12+K_s2
Dephosphorylation of Pjak	w_8=h4-x_13
Production of STAT	w_9=h5+x_14
Degradation of STAT	w_10=h5-x_14
Phosphorylation of STAT by pJAK	w_11=h6+x_14x_14+K_s3
Dephosphorylation of pSTAT	w_12=h6-x_15
Acetylation of pSTATn	w_13=h7+x_16
Deacetylation of pSTATn	w_14=h7-x_17
Inhibitory effect of pSTATnac on SHP1	w_15=h7i-x_17
Nuclear import of pSTAT	w_16=h_i8x_15
Nuclear export of pSTAT	w_17=h_e8x_16
Production of SOCS3 downstream the pathway	w_18=h9+x_16
Degradation of SOCS3	w_19=h9-x_18
Production of SMAD7 downstream the pathway	w_20=h10+x_16
Inhibitory effect of SOCS3 on STAT phosphorylation	w_22=h10i-x_18
Production of SHP1	w_23=h11+x_19
Degradation of SHP1	w_24=h11-x_19
Inhibitory effect of SHP1 on STAT phosphorylation	w_25=h11i-x_19
pSMAD3-pSTAT binding	z_1=g1+x_4x_15
pSMAD3-pSTAT unbinding	z_2=g1-x_20
Translocation of pSMAD3-pSTAT to nucleus	z_3=g_ix_20
Translocation of pSMAD3-pSTAT to cytoplasm	z_4=g_ex_21
Production of SNAIL in the pathway	z_5=g2+x_21
Degradation of SNAIL	z_6=g2-x_22

CXCL12 signaling pathway begins when this factor reaches its cell surface receptor, CXCR4. The binding of CXCL12 to its receptor prompts various gene expression and regulatory mechanisms through different intermediate proteins inside the cells. Phosphoinositide 3-kinases/ protein kinase B (PI3K/ PKB or Akt) mediator pathway initiate its activity in response to CXCL12 ligand-receptor phosphorylating action. It performs key functions such as controlling cell growth, proliferation, and migration by stimulating transcription factors such as Nuclear factor kappa-light-chain-enhancer of activated B cells (NFκB)[41,43-45]. Then, agitation of NFκB leads to overexpression of matrix metalloproteinases (MMPs) which are the major destructive factors of ECM[46]. ECM degrading or remodeling is a necessity for cancer cell invasion[47,48]. Thus, as it turns out, CXCL12-CXCR4 pathway is selected as a fundamental pathway for dynamical modeling of cancer cells. In our proposed model, CXCL12 pathway is governed by inducing CXCL12 and CXCR4 from CAFs and cancer cells, respectively.

LIF is the other abundant pro-invasive factor in the TME whose overexpression associates with poor clinical results for cancer patients. It displays strong upregulation as a target gene of TGFβ pathway. Moreover, it mediates fibroblasts reprograming into CAFs[36,37]. So, in the presented model, we consider LIF signaling pathway within CAFs in both paracrine and autocrine manners. Janus kinase/ signal transducer and activator of transcription (JAK/STAT) intermediate pathway is driven immediately after binding of LIF to its cell surface receptor. Then, the activated STAT enters nucleus and consequently contributes in transcription of target genes such as suppressor of cytokine signaling (SOCS) and SMAD7[37,49,50]. These two factors both regulate signal transduction by negative feedback creation in LIF and TGFβ pathways. As stated above, SMAD7 inhibits phosphorylation of SMAD2/3 and similarly STAT activation is prevented by SOCS protein[51]. There is another negative regulator for LIF pathway named small heterodimer partner (SHP) proteins which inhibits STAT phosphorylation[10]. Furthermore, if STAT protein undergoes the process of acetylation, it can reduce the expression level of SHP protein[10,37]. Recent studies have shown that SMAD3 and STAT3 cooperative positively to enhance epithelial–mesenchymal transition (EMT) via SNAIL upregulation[52]. EMT is an evolutionary process that confers invasion features upon tumor cells by changes in cellular phenotype[52,53]. It is worth noting that in LIF signaling cascade, SMAD7 upregulation and positive cooperation of SMAD3-STAT3 are crosstalk points between LIF and TGFβ pathways which are considered in our model. To the best of our knowledge, it is the first dynamic model which describes the latter crosstalk point in LIF and TGFβ pathways. Moreover, in order to make a comprehensive understanding the roles of TGFβ, LIF, and CXCL12 pathways in malignant tumors, the interactions mentioned above are mathematically illustrated in their dynamical model.

The dynamical systems unfolding the biochemical interactions in CAFs and cancer cells are based on nonlinear ODEs which characterize concentration changes of the involved molecules over time. Biochemical components of signaling pathways and their corresponding states in ODEs for cancer cells and CAFs are shown in Tables S1 and S2, respectively. The ODEs are extracted from laws of mass action and Michaelis–Menten which are used vastly to analyze reaction kinetics[54,55]. Mathematical representation of biochemical reactions involved in cancer cell and CAF dynamics are shown in Tables 1 and 2. Mathematical expression of a biochemical interaction includes multiple parameters as shown in these two tables. In Tables S3 and S4, model parameters, values and their references are provided. The values of kinetic parameters are gathered from literature whenever possible. Due to the large number of parameters and unavailability of their exact values in previous studies, most of the parameters in Tables S3 and S4 were estimated using HUKF approach. HUKF is an accurate and robust method designed for state and also parameter estimation of nonlinear systems that use discrete time data points as observed measurements. This algorithm relies on transformed statistics of parameters which are propagated through main steps including prediction and correction. In the first step, using model structure, the desired parameters are estimated and then priory estimates are optimized following the second step. Both model and observed data uncertainties are taken into account by the Kalman filtering approach which leads to reliable estimations[29,56]. In this study, we used two different gene expression profiling results datasets of time course studies (GSE17708 and GSE129189 for cancer cell and CAF, respectively) from Gene Expression Omnibus (GEO) database as observation data for HUKF algorithm[57,58]. Using these tables, cancer cell and CAF dynamics are modeled as state space systems of 21 and 26 nonlinear ODEs, respectively which can be found in Supplementary materials.

Mathematical representation of intercellular communications in the TME

The intercellular scale is defined by the communication signals exchanged between cancer cells and CAFs which allow for the description of cancer cells’ effecs on CAF functioning and consequently, the impact of CAFs on cancer cell invasion. The basic structure of the model presented in this section is inspired by our previous work in[29] with the notable exception that here we consider both spatial and temporal changes of signaling molecules in the model. As depicted in Fig. 4, TGFβ and LIF molecules transmit from cancer cells to CAFs where they stimulate the corresponding signaling pathways and lead to expression of metastatic factors such as CXCL12, SNAIL, TGFβ and LIF. Among these molecules, CXCL12 and TGFβ were identified as molecules that communicate signals from CAFs to cancer cells. They circulate in TME and once they reach the surface of cancer cells, tend to induce signal transduction.

Figure 4

Intercellular interactions between CAFs and cancer cells in the TME. TGFβ, LIF and, CXCL12 pathways are responsible for the release of communicating signals in the TME. TGFβ and LIF molecules are circulating from cancer cells to CAFs, as well as CXCL12 and TGFβ from CAFs to cancer cells.

To come up with a compendious view of signal transmission through the TME, we consider a three-dimensional space as tumor volume and its microenvironment. The mathematical representation for each signal is then provided using PDEs. Since gene expression of the signals, as well as their other regulatory activities occur within the cells, a number of terms in ODE systems of cancer cell and CAF will appear in PDEs.

On the other hand, when the transmission molecules circulating in the TME come into contact with the target cell, they are integrated onto the cell surface. The amount of reached factors on the cell is computationally modeled by an integral over the cell surface area. As a result, we have a coupled ODE–PDE framework which describes how CAFs and cancer cells interact in TME.

As it is shown in Table S1 and S2, our model includes four signaling molecules, . Spatial–temporal changes of the th molecule () is modeled using reaction–diffusion equation in (1):where and indicate first time derivative and Laplacian of diffusing signal ()‚ respectively. Also, is space vector and denotes the kinetics of interactions with other factors () in domain and time duration of where is the simulation step and is defined later. is a positive constant which represents diffusion coefficient of . The values of for each of the signals are shown in Table S5. In our model, is assumed to be a sphere with a radius of 2 cm[59,60] showing a three-dimensional space for the tumor. The Eq. (3) in Supplementary files shows the PDEs and complete form of dynamics for each signaling molecules, as well as their initial concentration. When solving a PDE, an initial condition and a boundary condition are typically needed, both of them are determined by the problem's time and space situations. We defined a constant value for initial condition () of each signal at the first step of simulation, [0,]. Then, at the end of each cycle, is updated and used for the next step. In addition, boundary condition on which shows the physical nature of signaling molecules on TME boundaries is considered as Dirichlet boundary condition.

As it is indicated in Fig. 4, signaling factors are integrated on each cell surface once they reached them and thus initiate the biological pathways within the cells. For determining the concentration level of molecules that is integrated () at the exact part of the domain (), we use the integral relation of (2):

In the next step of the intercellular level, (1) and (2) are used for describing the distribution of signaling molecules in TME. Then, with the exception of ODEs describing the signaling factors variations, we incorporate the obtained four PDE-integral systems for to with the previous dynamical systems. Eventually, a coupled ODE–PDE dynamical system emerges that is able to mathematically explain the conversations among a CAF and a cancer cell, as well as their intracellular dynamics at the intercellular level.

Biomechanical modeling movement of cancer cells in the TME

Cancer cells movement is occurred when they are prepared to detach from the tumor tissue and then leave the primary tumor under the influence of various intracellular and extracellular signal transduction pathways, as well as other TME components such as CAFs and ECM. ECM is one of the main TME components that is involved in metastasis. During the progression of metastasis, ECM is remodeled through a series of quantitative and qualitative modifications by degrading enzymes such as MMPs released by cancer cells[61]. So, in order to get a comprehensive view of cancer cells’ migration, we used an ODE adapted from[62] that included ECM generation and degradation to model the function of ECM in the TME as (3):where and are ECM and MMP concentrations, and and are generation and degradation coefficients of ECM, respectively.

We upgraded the model of intercellular interactions between one cancer cell and CAF to a population level as a necessary step before mathematical modeling of metastasis. Then, the status of each cancer cell was assessed by tracking each cell individually. In this case, each cancer cell is thought to be able to choose between quiescent, proliferative, or migratory states with regards to the concentration levels of certain effective biochemical factors including SNAIL, ECM and CXCL12. These factors are chosen because of their key roles in EMT, cancer cell adherence to the primary tumor, and tumor cells growth and proliferation, respectively[53,63,64]. The variations of CXCL12 and ECM in intercellular level of our model have already been modeled, but concentration of SNAIL in the TME must also be considered. Thus, we used (4) to mathematically represent how SNAIL is diffused in the TME following its upregulation in CAFs.where is the diffusion coefficient of SNAIL and its value is shown in Table S5.

As it is depicted in Fig. 5, the modeling framework for designating the status of cancer cells is based on the newly implemented coupled ODE–PDE dynamical system in the previous sections, as well as (3) and (4). During the time interval, the combined dynamical system is solved and the levels of ECM, SNAIL and CXCL12 factors are determined at each . Then, they are compared to a certain threshold value, and each cancer cell decides whether to keep or change its condition.

Figure 5

The multi-structure model for determining each cancer cell status. In this figure, u indicates the communicating signals which are transferred between cancer cells and CAFs in the tumor microenvironment (TGFβ and LIF molecules from cancer cells to CAFs, and CXCL12 and TGFβ from CAFs to cancer cells). Also, , , and represents CXCL12, ECM, and SNAIL concentration in the model. Furthermore, letter A indicates the surface area of a cancer cell or CAF which received a transduction signal and T implies the threshold amount of each signal which is indexed.

After assessing the status of cancer cells, we used an intuitive force-based biomechanical model based on the equation of motion to model movement of cancer cells. Cancer cells in migratory state will move in the population due to forces exerted in the TME as a result of interactions with other cells and the environment. In our model, cancer cells and CAFs are represented by spheres with radii of 8 and 7 µm, respectively[65,66] and it is assumed that the forces are applied to their centers. Also, the trajectory of each cell in the TME is specified by an equation of motion resulting in a Newtonian dynamic. The basic equation that governs the cancer cells movement is described as (5):where and are repulsive, and adhesive forces between each two cancer cells, respectively. Also, stands for haptotaxis force which is the directional motility of cells in response to gradients of adhesive surfaces such as ECM[67,68]. Moreover, denotes the active force that is driven by SNAIL concentration as the determining factor in metastasis. We also denote the friction between cancer cells, and the TME by [69]. Mathematical representations of these forces are shown in Table 3. More details about the constants and parameters of Table 3 can be found in Table S9. In each time interval , the velocity of cancer cells and consequently their displacement are calculated using (5) and Newtonian dynamic ().

Table 3

Mathematical representations of forces in the multiscale model.

Force	Mathematical representation	Reference
Repulsive force	Frep=43E^_iR^_iδij32	[69]
Adhesive force	Fadh=-πWR^	[69]
Haptotaxis force	F^hap=-χ∇ECM	[13]
Active Force	Fact=ηA^nA^n+A0nn^	[70]
Friction force	Ffric=-ΓcsV_	[18]

Experimental approaches

Cell preparation

MCF-7 cells were cultured in Dulbecco’s Modified Eagle Medium/Nutrient Mixture F-12 (DMEM/F12, Bio Idea; Iran) supplemented with 10% fetal bovine serum (FBS; Viva, Iran), 100 mg/ml streptomycin and 100 U/ml penicillin in 5% CO2 atmosphere at 37 °C. When cultured cells reached more than 90% of confluence degree, cell culture supernatants were collected and stored at − 80 °C.

To isolate breast cancer-associated fibroblasts (CAFs), breast tissues were minced in DMEM/F12 supplemented with 10% FBS, and 1% Penicillin/Streptomycin. The tissues were cut into small pieces (1–2 mm3) and cultured in the same media at 37 °C for about 2 weeks until Fibroblasts grew out of the tissue. Then, cells were removed and cultured in T25 flask at 37 °C in a humidified atmosphere with 5% CO2. When cancer-associated fibroblasts (CAFs) cells reached more than 90% of confluence, cell culture supernatant was collected and stored at − 80 °C.

Experiment

The CAFs were divided into two groups: group 1 (control, no treatment), and group 2 (cells treated with culture supernatant of MCF-7 cells). The CAFs were then incubated with DMEM F12 medium containing 10% FBS (250 µl) and supernatant collected from MCF-7 cells (750 µl) in a ratio of 1: 3 for 24 h, 48 h and 72 h. Also, MCF-7 cells were divided into two groups: group 1 (control, no treatment), group 2 (cells treated with culture supernatant of CAFs). MCF-7 cells were incubated with DMEM F12 medium containing 10% FBS and supernatant collected from CAFs in a ratio of 1: 3 for 24 h, 48 h and 72 h.

Wound scratch assay

Briefly, the MCF-7 cells, and CAFs (1 × 106 cells/well) were plated in 12-well plates for 48 h to a confluence of about 90%, then wounded by scratching with a p100 pipette tip. Then, the debris was removed and the cells were washed once with 1 mL of the growth medium to assure the edges of the scratch were smoothed by washing. As mentioned above, the cells were cultured for 24 h, 48 h and 72 h. After different time courses, cells were washed twice with PBS and the wound was observed under a microscope. Cells were photographed before and after 24, 48 and 72 h of incubation and images were analyzed with ImageJ software to calculate the area of each scratch. Images of each well at different times were compared to comprehend the effect of treatments on cell migration[71].

RNA extraction and cDNA synthesis

As mentioned above, the cells were cultured for 24 h, 48 h and 72 h. After each time course, total RNA was extracted from CAFs, and MCF-7 cells using the BioFACT™ Total RNA Prep Kit (BioFACT, South Korea), according to the manufacturer’s recommendations for cultured cells. RNA quantity was assessed using NanoDrop spectrophotometer. The average RNA yield was 0.5 μg RNA per 106 cells. RNA was subsequently reverse-transcribed to cDNA using the BioFact™ RT Series (BioFACT, South Korea) according to the manufacturer’s recommendations. Then, they were incubated for 60 min at 42 °C and were terminated the reaction for 5 min at 95 °C.

Real-time polymerase chain reaction (PCR)

Quantitative PCR was performed in a Rotor-Gene 6000 (Corbett Life Science, Sydney, Australia). For mRNA quantification, the 2X Real-Time PCR Master Mix (BioFACT™, High ROX, Korea) was used in combination with pre-designed primers for Leukemia inhibitory factor (LIF), Transforming growth factor beta (TGF-β), SMAD7, C-X-C motif chemokine 12 (CXCL12) and for the reference gene GAPDH as the internal control, were designed using the Allele ID software and BLAST (NCBI online server).

All real-time PCR reactions were performed in duplicate at a final volume of 10 μl per well using a qPCR Master mix (Jena Bioscience, GmbH) and a Rotor-Gene 6000 (Corbett Life Science, Sydney, Australia). The following thermal cycling conditions were applied: polymerase activation/denaturation at 95 °C for 15 min, and 45 amplification cycles at 95 °C for 20 s, 60 °C for 20 s and 72 °C for 30 s. The mRNA fold increase or fold decrease with respect to control well was determined and Relative quantification (ΔΔCT) method was employed. The relative quantification of mRNAs was calculated using the 2 −ΔΔct method and according to the formula normalization ratio (N.R.) = 2−ΔΔct.

Ethics declaration

We confirmed that all experiments were performed in accordance with relevant guidelines and regulations. Cancer-associated fibroblasts (CAFs) were isolated from surgical specimens taken from the tissue of patients with breast cancer at Seyed-al-Shohada Hospital, Isfahan. A written informed consent was taken from all participants. This study was approved by Isfahan University of Medical Sciences (the registration number: IR.MUI.MED.REC.1398.727).

Results

Nonlinear ODE modeling of cancer cell and CAF dynamics explains SMAD7, TGFβ, LIF and CXCL12 gene expression profiles

Mathematical modeling of intracellular interactions in cancer cells and CAFs allows determining gene expression variations of biochemical species involved in the reactions. The models presented in supplementary files were simulated using MATLAB ODE toolbox 2020a[72]. Figures 6 and 7 indicate the expression of fate-determining factors such as TGFβ, LIF, and CXCL12 increases in both cancer cells and CAF throughout the first three time periods. Also, simulation results succeeded to reproduce the behavior of SMAD7, TGFβ, LIF, and CXCL12 molecules in both cancer cell and CAF. This may be explained by looking at the goodness of fit criteria in Tables S5 and S6, which show how well our model outputs explain the variation in experimental data. Moreover, Bland–Altman graphs were provided for output genes as supplementary figures (Figs. S1 and S2) to highlight the agreement between two different datasets (model output and experimental data). Generally, the reported values and graphs for performance evaluation of output genes over four time points indicate a close enough proximity for model validation. Furthermore, using Morris method, we performed a global sensitivity analysis to evaluate the parameter effects on the model outputs[73]. Due to the large number of parameters, we selected a set of parameters which are directly affect the model outputs. In this method, selected parameters are changed over the entire range of probable values and the differences of the desired elements during these variations are calculated. Then for each parameter, mean and variance values are computed which are the indicators of effectiveness and nonlinearity of the parameter on the outputs, respectively. Sensitivity analysis showed that expression of TGFβ and its receptor, as well as expression of SMAD7 in this pathway are among the principle reactions that have the biggest effect on the outputs. More details on the selected parameters and also their mean and variance values obtained during the sensitivity analysis could be found in the Supplementary files.

Figure 6

Experimental data and model outputs of ODE modeling of cancer cell. Gene expression profiles of SMAD7, TGFβ, LIF, and CXCL12 molecules in MCF-7 cell over four time points (0, 24, 48, and 72 h) and their corresponding behavior in the dynamic model are depicted by red dots and black curves, respectively.

Figure 7

Experimental data and model outputs of ODE modeling of CAF. Gene expression profiles of SMAD7, TGFβ, LIF, and CXCL12 molecules in CAF over four time points (0, 24, 48, and 72 h) and their corresponding behavior in the dynamic model are depicted by red dots and black curves, respectively.

A force-based scheme incorporating the coupled ODE–PDE model and biomechanical forces describes cancer cells’ movement in population level

In order to investigate the cellular level of communication between cancer cells and CAFs, several PDEs are combined with ODE dynamics of each individual cell to model the interacting signals between cancer cells and CAFs circulating within the TME. The coupled ODE–PDE system is solved using finite element method (FEM) in MATLAB PDE toolbox 2020a. FEM is a numerical technique to solve PDEs over a defined space by dividing it into smaller parts[74]. The primary tumor in our model is a sphere with a radius of 2 cm that has been preprocessed by meshing to around 1000 subdomains. Also, the position of cancer cells and CAFs are randomly determined in the tumor. Then, using PDEs described as (1), TGFβ and LIF signals from cancer cells as well as TGFβ and CXCL12 molecules from CAFs release into the TME. Immediately after reaching these signals to the opposite cell type, signaling pathways within each cell get started according to the dynamics of cancer cell and CAF during a time interval. Furthermore, initial values of the next time step are updated in each . So, the cellular conversations between cancer cells and CAFs in a three dimensional domain are modeled using a coupled ODE–PDE dynamic system.

Then at the population level, the model has been developed to display movement of cancer cells as a result of their interactions with TME components and the designated status of each cell. In this scenario, cancer cells in migratory condition are under the influence of biomechanical forces that result from local interactions with neighbor cancer cells, CAFs and cancer cells. In summary, after solving coupled ODE–PDE dynamical system in a time interval, the level of diffused signals on each cell at are identified. Then, status of each cancer cell is determined in one of the quiescent, proliferating or migratory condition by the instruction illustrated in Fig. 5. Eventually, based on (5) and the defined forces in Table 3, displacements of cancer cells in migratory status are calculated. These events are repeated every seconds and at the end of each time step, initial values are updated for the next modeling step.

Starting with 20 CAFs and 20 cancer cells randomly placed in a sphere as the primary tumor, this multiscale model is simulated using MATLAB 2020a. The process is repeated every and the results are depicted in Fig. 8. Initial cancer cells in the TME are considered to be in quiescent state and their primary arrangement is randomly selected. To validate the model, a wound scratch assay has been performed to evaluate MCF-7 cells migration while they are co-cultured with CAFs and the movement are assessed using Image J software. Comparison of the average displacement in wound scratch assay and model outputs are shown in Table 4, and the performance evaluation criteria in Table S8, as well as Bland–Altman graph in Fig. S3, indicate a proper performance of the applied techniques to model cancer cell metastasis.

Figure 8

Simulation results of multiscale model of cancer cells movement during 10 time steps. Cancer cells and CAFs are depicted by red and black dots, respectively. Primary arrangement of cancer cells is randomly selected and it is shown in k = 0. After 10 repetitions, not only has the number of cells increased, but also they have been moved in space and left the primary tumor. The x–y-z axes have been scaled in centimeters.

Table 4

Average displacement of cancer cells in wound scratch assay versus model outputs.

Time (hour)	Displacement (mm) in wound scratch assay	Displacement (mm) in the model
0	0	0
24	0.158	0.180
48	0.223	0.264
72	0.672	0.634

Discussion

It is well understood that cancer cells’ invasive features are highly dependent on their interactions with TME components, and thorough understanding of the communications could contribute to the design of new therapeutic approaches for preventing metastatic progression. Biological pathways, transmission signals, and mechanical forces occurring at different physical scales all influence the migration of cancer cells. In this paper, we have presented a multiscale model of metastasis with an emphasis on the impact of TME components such as CAFs and ECM. The model consists of three principle regulatory scales including intracellular, intercellular, and population levels which are validated against experimental data.

In the first step, biological pathways within the cancer cells and CAFs are modeled using ODE-based dynamical systems. TGFβ and CXCL12 pathways, as well as their crosstalk have been considered as the main regulatory mechanisms of cancer cell dynamics. In addition, an ODE model of CAF dynamic is derived using the mathematical depiction of physiological events involved in TGFβ and LIF signaling pathways. Since there are no exact values for several number of kinetic parameters, they are estimated using the HUKF algorithm. Moreover, the dynamical models produced results that are very similar to the experimental data of SMAD7, CXCL12, LIF, and TGFβ gene expression profiles. In this study, two different sets of experimental data have been used for the procedure of parameter estimation and model validation. Despite the fact that model yielded results in a close agreement to experimental data, small number of data points could be considered as a limitation of the study. However, it is inevitable that mathematical modeling of biological phenomena might be compelled to apply the limited number of data points, due to this fact that determining a large number of experimental data points can be impractical or impossible for a variety of reasons. To ensure meaningful results, we used technical methods such as HUKF algorithm which considers the uncertainty not only in the model parameters but also in the measurement datasets.

In order to describe how LIF, CXCL12, and TGFβ molecules which are considered as communicating signals circulate within the TME, we have considered a system of PDEs in a three-dimensional domain. The PDE system depicts the concentration variations of signaling molecules over time and space. Then, ODE models of cancer cells and CAFs are coupled to these PDEs through the mathematical statements utilized for evaluating the amount of signals received on the cell surfaces. As a result, a coupled system will be developed that can demonstrate the conversations of CAFs with cancer cells as well as their effect on tumor progression and invasion. Furthermore, by combination of coupled ODE–PDE modeling framework with ECM characteristics, the model succeeded to determine the status of each cancer cells in one of quiescent, proliferating or migratory conditions. Consequently, proliferating and migratory cells begin to replicate and migrate, causing tumor growth and invasion, respectively. Cancer cells’ movement are regulated by an equation of motion that includes adhesive, repulsive, haptotaxis, active, and friction forces applied to cancer cells. To demonstrate the functionality of the multiscale model, we performed a wound scratch assay test for measuring the displacement of MCF-7 cells co-cultured with CAFs at specific time points, and favorably acceptable accordance is observed when measured displacement in wound scratch assay and model outputs are compared. Although there are various methods of measuring directional collective cell migration, this assay was performed for measuring the tumor cells displacement due to its versatile application to fibroblasts and tumor cell types, simple and fast protocol, and easy analysis method. Moreover, not also cell migration, but also cell proliferation contribute to the wound healing, and these phenomena both has been considered in our modeling framework[96].

Despite the fact that a number of experimental studies looked into the function of CAFs in malignant tumor metastasis, the mechanisms behind their influences on the invasiveness of cancer cells remained elusive, especially in the case of their mutual interactions in the TME. The proposed mathematical model is able to provide a comprehensive framework to illustrate not only the internal dynamics of CAFs and cancer cells, but also their interactive conversations throughout the TME. Furthermore, regulatory feedback loops that play a vital role in cancer cell and CAF phenotype change have been considered in the current study, which have not been investigated quantitatively so far. The model also shed more light on motility of cancer cells by taking into account ECM remodeling and its critical features in EMT. ECM as a TME component provides structural support to surrounding cells, but during advanced stages of metastasis, it undergoes some modifications and fails to keep cancer cells adhered together. As a result, cancer cells detach from the tumor tissue and leave the primary tumor. Moreover, the ability of the model to designate the cancer cells states makes an alive structure that is proliferating and moving, as well as allowing us to track the behavior of each cell during the simulation. Therefore, considering the fact that findings of the current multiscale model have been confirmed by experimental evidence, this paper can be regarded as a pioneer in the study of complicated phenomena such as tumor metastasis. From both empirical and computational point of view, our approaches pave the way for further studies and designing more therapeutic methods in cancer cell metastasis.

Supplementary Information.