Research article Special Issues

Mathematical modeling and numerical simulation of a multiscale cancer invasion of host tissue

  • In this study, we develop an advection-reaction-diffusion system of partial differential equations (PDEs) to describe interactions between tumor cells and extracellular matrix (ECM) at the macroscopic level. At the subcellular level, we model the interaction of proteolytic enzymes and the ECM with a set of ordinary differential equations (ODEs). A contractivity function is used to couple the macroscopic and microscopic events. The model is supplemented with nutrients supply from the underlying tissue. These PDE-ODE systems of equations model the on-set of tumor cell invasion of the host extracellular matrix. The model accounts for different time and spatial scales at the macroscopic and microscopic levels. Contact inhibition between the tumor cells and the tumor micro-environment are also accounted for through a nonlinear density-dependent diffusion and haptotaxis coefficients. In the numerical simulations, we use a nonstandard finite difference method to illustrate the model predictions. Qualitatively, our results confirm the three distinct layers of proliferating, quiescent and necrotic cells as observed in multicellular spheroids experiments.

    Citation: Peter Romeo Nyarko, Martin Anokye. Mathematical modeling and numerical simulation of a multiscale cancer invasion of host tissue[J]. AIMS Mathematics, 2020, 5(4): 3111-3124. doi: 10.3934/math.2020200

    Related Papers:

    [1] Ali Sadiq Alabdrabalnabi, Ishtiaq Ali . Stability analysis and simulations of tumor growth model based on system of reaction-diffusion equation in two-dimensions. AIMS Mathematics, 2024, 9(5): 11560-11579. doi: 10.3934/math.2024567
    [2] Juya Cui, Ben Gao . Symmetry analysis of an acid-mediated cancer invasion model. AIMS Mathematics, 2022, 7(9): 16949-16961. doi: 10.3934/math.2022930
    [3] Cagnur Corekli . The SIPG method of Dirichlet boundary optimal control problems with weakly imposed boundary conditions. AIMS Mathematics, 2022, 7(4): 6711-6742. doi: 10.3934/math.2022375
    [4] Fawaz K. Alalhareth, Seham M. Al-Mekhlafi, Ahmed Boudaoui, Noura Laksaci, Mohammed H. Alharbi . Numerical treatment for a novel crossover mathematical model of the COVID-19 epidemic. AIMS Mathematics, 2024, 9(3): 5376-5393. doi: 10.3934/math.2024259
    [5] Hsiu-Chuan Wei . Mathematical modeling of ER-positive breast cancer treatment with AZD9496 and palbociclib. AIMS Mathematics, 2020, 5(4): 3446-3455. doi: 10.3934/math.2020223
    [6] Yumei Chen, Jiajie Zhang, Chao Pan . Numerical approximation of a variable-order time fractional advection-reaction-diffusion model via shifted Gegenbauer polynomials. AIMS Mathematics, 2022, 7(8): 15612-15632. doi: 10.3934/math.2022855
    [7] Kaihong Zhao . Attractor of a nonlinear hybrid reaction–diffusion model of neuroendocrine transdifferentiation of human prostate cancer cells with time-lags. AIMS Mathematics, 2023, 8(6): 14426-14448. doi: 10.3934/math.2023737
    [8] José L. Díaz . Existence, uniqueness and travelling waves to model an invasive specie interaction with heterogeneous reaction and non-linear diffusion. AIMS Mathematics, 2022, 7(4): 5768-5789. doi: 10.3934/math.2022319
    [9] Anusmita Das, Kaushik Dehingia, Hemanta Kumar Sharmah, Choonkil Park, Jung Rye Lee, Khadijeh Sadri, Kamyar Hosseini, Soheil Salahshour . Optimal control of effector-tumor-normal cells dynamics in presence of adoptive immunotherapy. AIMS Mathematics, 2021, 6(9): 9813-9834. doi: 10.3934/math.2021570
    [10] Iqra Batool, Naim Bajcinca . Stability analysis of a multiscale model including cell-cycle dynamics and populations of quiescent and proliferating cells. AIMS Mathematics, 2023, 8(5): 12342-12372. doi: 10.3934/math.2023621
  • In this study, we develop an advection-reaction-diffusion system of partial differential equations (PDEs) to describe interactions between tumor cells and extracellular matrix (ECM) at the macroscopic level. At the subcellular level, we model the interaction of proteolytic enzymes and the ECM with a set of ordinary differential equations (ODEs). A contractivity function is used to couple the macroscopic and microscopic events. The model is supplemented with nutrients supply from the underlying tissue. These PDE-ODE systems of equations model the on-set of tumor cell invasion of the host extracellular matrix. The model accounts for different time and spatial scales at the macroscopic and microscopic levels. Contact inhibition between the tumor cells and the tumor micro-environment are also accounted for through a nonlinear density-dependent diffusion and haptotaxis coefficients. In the numerical simulations, we use a nonstandard finite difference method to illustrate the model predictions. Qualitatively, our results confirm the three distinct layers of proliferating, quiescent and necrotic cells as observed in multicellular spheroids experiments.


    Mathematical modeling of biological systems especially on cancer has in recent times received much attention because, mathematical models of relevant biological processes and its associated numerical simulation, has reduced the complicated and costly experimental procedures.This is buttressed by the fact that simulation of cancer behavior across multiple biological scales in space and time is increasingly being recognized as a powerful tool that enables more accurate predictions, see [1,2,3,4,6,23] and references therein. The primary objective of multidisciplinary efforts in cancer research from the fields of mathematics, physics, and chemistry just to mention a few, is to understand the etiology of the disease and to develop suitable diagnosis and treatment.

    Cancer is a generic term given to a group of diseases in which abnormal cells grow uncontrollably. Additionally, cancer cells can invade surrounding tissues and further spread to other organs in the body (http//www.who.int/health-topics). Cells become cancerous from a series of molecular (gene mutations) events that alter their fundamental properties and structural rearrangements. The resultant effects of this structural rearrangement is a characteristic abnormality in the process regulating the cells carefully regulated order of proliferation, differentiation, and death. If this abnormality continues after many years, the cells evolve through different premalignant stages to become invasive and degrade the extracellular matrix (ECM) or surrounding tissue. Such abnormal cells proliferate exponentially to form a lump called a tumor. There are two types of tumors: benign and malignant tumors, and whiles the former is not harmful the later is cancerous. Malignant tumor cells grow uncontrollably and may spread to other organs and neighboring tissues in the body through the metastasis process. Cancer development has three distinct stages: avascular, vascular, and metastatic [22]. A review of the literature shows that though all the three stages have been well studied and documented, work on avascular tumors dominates the number of research among the three. This does not in any way imply the importance of avascular tumors in cancer research but rather it is the simplest to model mathematically and to validate using experimental results (eg. from multicellular spheroids). Secondly, the worth of information drawn from research on avascular tumors can increase our understanding of the last two stages (vascular and metastasis). It is worth mentioning that the invasion begins as soon as the tumor cells can penetrate the surrounding tissue. Cells in the avascular tumor stage are not cancerous but form a lump due to some abnormality in the growth and agammaegation of proliferating cells.

    Though the early stages of tumors are difficult to study due to their small size, education on the disease and technological advancements have contributed immensely to the study of the early stages of the tumor through in vitro experiments using multicellular spheroids, -an experimental process of growing seed cell(s) taken from a tumor cell line in a liquid medium supplemented with sufficient nutrients, see [13,20,21].

    Tumor cells receive nutrients (eg., oxygen and glucose) supply from surrounding tissues. These nutrients reach the central core of the tumor by diffusing through layers of cells. Therefore, as the tumor increases in size, the amount of nutrients from the tumor microenvironment (TME) reaching cells at the inner rim (center) decreases from the outer layer until it becomes insufficient to sustain viable cells at the center of the tumor. Experimental studies of multicellular spheroids exhibit three distinct characteristic layers of nutrient-rich proliferating cells and an inner layer of dead (necrotic) nutrient-starved cells. These two layers are separated by a layer of nutrient deficient quiescent (hypoxic) cells. Quiescent cells are alive and can diffuse but are non-proliferating. Cells in the quiescent state may proliferate when they receive sufficient nutrients. The tumor size may be limited by the amount of nutrients that reach the inner layer of the tumor. Thus as the tumor grows in size, the volume of necrotic cells at the central core also increases. Tumor growth stops when the balance of cell proliferation and cell death is reached [15,22]. Therapies for cancer treatment usually target proliferating cells therefore, the population of quiescent cells may be liable to treatment failures [20,21,22]. It is therefore very crucial to explore and investigate the population of quiescent cells as we report in this study.

    Insufficient supply of nutrients into the inner rim of the lump of cells is a contributing factor that causes the cells to initiate a process to acquire their blood supply (angiogenesis). Once the cells succeed to acquire their blood supply, the tumor cells can escape from one primary site to a secondary site to initiate secondary tumors through the circulatory system (metastasis) at other parts of the body. At the angiogenesis stage, the tumor is vascular and there is a decrease in the cell-cell and cell-matrix adhesion [1]. Vascular and metastatic stages are those that reduce the quality of an individual patient's life.

    In furtherance to cell-cell and cell-ECM adhesion, cell surface receptors bind various ECM fibers to enable cell-matrix adhesion [12]. By adhering to the ECM fibres, cells can migrate from a region of lower concentration of ECM adhesive fibers to a region of higher concentration of ECM adhesive fibers in a process known as haptotaxis. Furthermore, cells migrate in response to nutrients/chemoattractant (or chemorepellent) in a process known as chemotaxis. See a detailed discussion of these processes in [3,6,10].

    The past five years have seen tremendous in vitro and in vivo research efforts aimed at understanding the etiology of cancer cell's tissue invasion and to develop appropriate therapies for treatment. Mathematical models on cancer have either been discrete or continuum for cell dynamics [33]. The first mathematical model for tumors formulated in terms of densities of live and dead cells was proposed by Ward and Kind (1997, 1999), in contrast to earlier models which were formulated based on discrete cell populations. Though clinically, the latter formulation is artificial and not supported by experimental evidence. In Ward and King (1997), tumor growth was entirely due to the generation of new cell division. The analysis of interdependence between viable proliferating cell rim, growth velocities, and model parameters were also reported in their model.

    In [11], the authors report that the relationship between solid tumors, ECM and multiple cell types in their TME enables the solid tumor to transform systematically into an invasive, metastatic and lethal disease.

    Cancer cells secrete proteolytic enzymes (metalloproteinases (MMP) and or urokinase plasmogen activator (uPA)) which bind the ECM fibers to transform their conformation and configuration. Furthermore, the binding can break the natural and physical barriers that restrict the cancer cells from spreading into surrounding tissues. Mathematical models for the biologically multiscale character of cancer invasion of host tissues comes with challenges emanating from the complex interplay between the tissue-scale (macro) and the subcellular, microscale events involving the proteolytic enzymes binding to the ECM, see [11,18,24,25,26,27]. The micro and macro scales differ in space and time as we discuss in section 2. The ECM acts as a scaffold for cells as well as a signaling pathway through which the cells communicate specific intracellular and extracellular instructions for cell proliferation, migration, differentiation and most importantly, the secretion of proteolytic enzymes that destroy the ECM [31,32]. In their two-scale moving boundary model, the authors in [11] report on the effect of proteolytic enzymes (uPA) at the molecular level on the tissue-scale (macro) dynamics at the edge of the tumor invasion. Further research on cancer invasion of host tissues and investigation on effect of uPA and MMP on the evolution of cancer cells through mathematical modeling backed by experimental data include [3,11,28,29,30]. Chaplain and Lolas [28] used a system of reaction-diffusion-taxis (RDT) partial differential equations to explore and report on the role of uPA's in the invasion, growth and spread of cancer cells whiles Chaplain and Deakins [30] studied the role of MMP's in cancer growth and spread. The proof of global existence and unique classical solution for the system of RDT PDE's using a priori estimates to the model proposed by Chaplain and Lolas [6] is given in [16]. Further models for cancer invasion include [1,2,3,4,5,6,7,8,9,10,17], and references therein.

    The paper is structured as follows: In Section 2 we discuss the formulation of the mathematical model, followed in Section 3 by a nondimentionalization and scaling of the model. In Section 4, detailed discretization of the model, special treatment of the nonlinear diffusion and haptotaxis terms and numerical solutions of the model at specific time intervals is presented. In Sections 5 and 6 we discuss the findings and conclusion of the model respectively. }

    The study focuses on the first step of cancer invasion of ECM. We use a continuum approach to capture the entire micro and macro scale volumetric tumor growth dynamics. The model has two parts: first, a set of ordinary differential equations that models the biochemical dynamics of the tumor cells interaction with ECM at the subcellular (micro) level. In the second part, a system of partial differential equations is used to model the evolution of proliferating cells, quiescent cells, necrotic cells, and ECM densities at the macroscopic level. Unless otherwise mentioned the terms tumor cells and cancer cells are used interchangeably. The above can best be described as an advection-reaction-diffusion (ARD) system of equations. In the next subsection, we present the events at the microscopic level followed by the macroscopic level events.

    Cancer cell invasion of tissues involves biochemical processes taking place at the subcellular, microscopic-scale level. These subcellular level events affect the cell motility and haptotaxis of tumor cells and intrinsically the invasion of ECM. The events at the subcellular level take place on scales in space and time that are different from the macroscale events. To observe the effects of the subcellular events at the macroscopic level, there is a time lag between the biochemical event at the subcellular level and when it is observed at the macroscopic (tissue) level. We will account for the time lag in a contractivity function used to couple the subcellular, microscopic scale events to the macroscopic scale events. It is important to mention that in between these two time-scales (micro and macro) is an intermediate time-scale known as the meso-scale (not discussed in this work), see [3] for more details. From the law of mass action, we can summarize the reaction between the proteolytic enzymes (l) and the ECM fibers (v) in a reversible reaction mechanism in Eq. (2.1). Equation (2.1) holds since none of the reacting proteins and products formed is volatile, precipitate or is consumed in other reactions.

    l+vk+klv (2.1)

    where k+ and k are association and dissociation constants respectively. Let l0 and y, denote the initial concentration of proteolytic enzymes and the concentration of proteolytic enzymes bound to ECM proteins (lv) respectively, then it holds l0R+. The concentration of unbound proteolytic enzymes on the cell surface is given by l0y. An ordinary differential equation (ODE) for the bound proteolytic enzymes is given in Eq. (2.2).

    dydθ=k+(l0y)vky (2.2)

    where θ is the microscopic time scale.

    It is assumed that the onset of tumor invasion of host tissue is predominantly dominated by random motion, and haptotaxis due to spatial gradients as a result of ECM fiber degradation. At the macroscopic level, we formulate our model in terms of a continuum of densities of necrotic cells, quiescent cells, proliferating cells and ECM fibers denoted by u(x,t), q(x,t), p(x,t) and v(x,t), respectively. The model is supplemented with nutrient concentration c(x,t) from the underlying tissue. The spatial coordinate in the one-dimensional domain is x and t is time. By this formulation, our model is geared towards an in vivo than an in vitro as this approach also allows us to account for the onset of invaded tissue in the model.

    As in [6,34], a conservation equation for these two contributions to the tumor cell flux Jn is given by Jn=Jrandom+Jhapto.. We assume that proliferating and quiescent cells migrate whiles necrotic cells are immotile. The diffusion of cancer cells into the surrounding tissues is given by Φj(κ,p,q,v)j, j{p,q}, where the non-linear density-dependent diffusion function Φj(κ,p,q,v) depends on the contractivity function κ, proliferating and quiescent cell densities (accounting for contact inhibition of migrating cells) and the ECM density v(x,t). The cells need to adhere to adhesion sites located on the ECM fibers for migration. It is assumed that the diffusion of proliferating and quiescent cell population are equal. Therefore, we use a non-linear diffusion function Φq(κ,p,q,v)q for quiescent cells, similar to the diffusion term describing proliferating cells [14,15]. The degradation of ECM fibers by proteolytic enzymes creates spatial gradients which directs the migrating cancer cells. Therefore, the function Ψ(κ,v)pv describes the haptotactic movement of cancer cells when they bind to the adhesion sites on the ECM [16]. Proliferating cancer cells compete for locally available space with the ECM fibers. This is described by the last but one term in the first equation (2.3) with a proliferating rate g(c) that depends on the concentration c(x,t) of available nutrients. When proliferating cells become deficient in nutrients, they move into a quiescent state. This is modeled by the term f(c)p. f(c) describes the rate at which proliferating cells transition into the quiescent state. The function h(c) describes the rate at which quiescent cells die due to nutrients starvation. See the explicit form of the functions f(c) and h(c) as given in [15], below. In the model, ECM is non-diffusive, we assume that proliferating and quiescent cells come into contact with the ECM to degrade it at a rate v. The ECM fibers reestablish themselves at a rate μv, whiles competing for locally available space with the cells. The competition is modeled by a term similar to that of the proliferating cells. Necrotic cells do not diffuse, rather they are produced when quiescent cells die. A contractivity function κ is used to couple the subcellular, microscopic scale biochemical events to the macroscopic scale (population of cells) events. This function determines the strength of tumor cell motility. The contractivity function has a constant time lag (τ) to account for the delay in the transfer of signals from the subcellular events to affect the dynamics of motile tumor cells on the macroscale. The complete system of equations describing the interactions between proliferating cells, quiescent cells, necrotic cells, ECM, proteolytic enzymes, contractivity function and cell nutrients are given by:

    {tp=(Φp(κ,p,q,v)p)(Ψ(κ,v)pv)+g(c)p(1r3vKv)f(c)p,t>0,Ωtq=(Φq(κ,p,q,v)q)+f(c)ph(c)q,t>0,Ωtu=δuq,t>0,Ωtv=δv(p+q)v+μvv(1r3vKv),t>0,Ωtc=Dc2c+k1c0(1αKr(p+q))k1ck2pc,t>0,Ωθy=k+(l0y)vky,θ>0,Ωθκ=aκ+My(tτ),θ>0,Ω. (2.3)

    where tc=ct for the PDEs, ΩRd for d=1 representing a bounded domain with a smooth boundary. We supplement the model in (2.3) with boundary and initial conditions, in sections 2.3 and 2.4. Whereas g(c) is an increasing function, f(c) and h(c) are both decreasing functions. The rate at which proliferating cells become quiescent is greater than the rate at which quiescent cells die to become necrotic and therefore, f(c)>h(c). Kr and Kv are carrying capacities of cancer cells and ECM respectively, δu is the rate at which cells become necrotic, α, k1 and k2 are constants, and τ>0 is a fixed constant time delay. Additionally, g(c)=βexpβc, f(c)=0.5(1tanh(4c2)) and h(c)=f(c)2. The non-linear density-dependent diffusion and haptotaxis functions are given by:

    Φp(κ,p,q,v)=Dpκ1+(pp+q)vKv, Φq(κ,p,q,v)=Dqκ1+(qp+q)vKv,

    Ψ(κ,v)=DhκvKv+v, r3=(pKr+qKr+uKr),

    p.υ=q.υ=u.υ=v.υ=c.υ=0 on (0,T)×Ω (2.4)

    where υ denote the outward unit normal vector on Ω.

    p(0,x)=p0(x), q(0,x)=q0(x), u(0,x)=0, κ(0,x)=κ0(x),

    c(0,x)=c0(x), y(0,x)=y0(x), t(,0], xΩ

    We rescale the system of equations in (2.3) using p:=pKr, q:=qKr, v:=vKv, u:=uKr, c:=cC0, x:=xL, t:=tT, t:=ηθ and y:=yl0. L and T are reference length and time scales respectively. We obtain the nondimensionalized system of Eq. (3.1), where the asterix has been dropped for convenience.

    {tp=.(Dpκ1+pp+qv).(Dhκv1+vpv)+λ(g(c)p(1(p+q+u)v)f(c)p),t>0,Ωtq=.(Dqκ1+qp+qv)+λ[f(c)ph(c)q],t>0,Ωtv=δv(p+q)v+μvv(1r3v),t>0,Ωtu=δuq,t>0,Ωtc=DcΔc+k1(1α(p+q))k1ck2pc,t>0,Ωθy=k+(1y)vky,θ>0,Ωθκ=aκ+My(θτ),θ>0,Ω. (3.1)

    where

    Dp=DpTL2, Dh=DhTKvL2, Dq=DqTL2, λ=1T.

    δv=δvTKr, μv=μvT,

    δu=δuT

    Dc=DcTC0L2, k1=k1T, k2=k2KrT

    k+=k+ηTKv, k=kηT

    a=aηT, M=MηT, τ=τηT.

    The equations are discretized using the implicit-explicit finite difference method. The space interval [0,1] is divided into k parts giving us k+1 nodes. The space step Δx=0.01. We use forward difference for the time derivatives and central difference for the diffusion and haptotaxis terms. Eg. ut=un+1unΔt, ux=ui+1ui12Δx and 2ux2=ui+12ui+ui1(Δx)2. To simplify the model for discretization purposes we define the following variables: Sp=pp+q, Sq=qp+q, M=κ1+Sjv and N=κv1+vp, where j{p,q}. The evolution equation describing the dynamics of the proliferating cell from equation (3.1) then reads:

    pt=.(DpMp).(DhNv)+λg(c)p(1(p+q+u)v)λf(c)p. (4.1)

    Systematically, the system of equations are discretized and numerically solved. It is important to note that the micro-level events update the contractivity function. As mentioned earlier, there is no spatial movement in the equation for ECM density therefore, the ODE is discretized using forward difference in time as follows:

    vn+1ivniDt=δv(pni+qni)vn+1i+μvvni(1rnivni)
    vn+1i=11+Dtδv(pni+qni)(vni+Dtμvvni(1rnivni)),i=0,1,2,3,...k

    where n denote the time level, Δt=0.01 is the time step for the microscale events and Dt=ηΔt is the time step for the macroscale events with η as the time scaling factor. The equation describing proteolytic enzyme interactions with ECM and the contractivity functions are also discretized respectively as:

    yn+1i=11+Δtk+vn+1i+Δtk(yni+Δtk1vn+1i)κn+1i=11+aΔt(κni+(ΔtMyn+1i))

    where k+=2, k=0.6 and a=3. In discretizing the PDE for cancer cells, we use the nonstandard finite difference scheme employed in [3]. By this approach, the diffusion and haptotaxis coefficients of the equation for the proliferating cells are handled explicitly while the remaining terms are handled implicitly. The discretization is given below

    pn+1ipniDt=Dp[Mp+MΔp]Dh[Nv+NΔv]+λg(c)p(1(p+q+u)v)λf(c)p=Dp(Mn+1i+1Mn+1i1)2Δx(pn+1i+1pn+1i)Δx+DpMni(pn+1i+1pn+1i)+(pn+1i1pn+1i)(Δx)2Dh(Nn+1i+1Nn+1i1)(vn+1i+1vn+1i)2(Δx)2DhNni(vn+1i+1vn+1i)+(vn+1i1vn+1i)(Δx)2+λg(cni)pni(1(rni)vni)λf(cni)pn+1i=Dp2(Δx)2(Mn+1i+1+Mni)(pn+1i+1pn+1i)+Dp2(Δx)2(Mn+1i1+Mni)(pn+1i1pn+1i)Dh2(Δx)2(Nn+1i+1+Nni)(vn+1i+1vn+1i)Dh2(Δx)2(Nn+1i1+Nni)(vn+1i1vn+1i)+λg(cni)pni(1(rni)vni)λf(cni)pn+1i=Dp2(Δx)2kRi(Mn+1k+Mn+1i)(pn+1kpn+1i)Dh2(Δx)2kRi(Nn+1k+Nn+1i)(vn+1kvn+1i)+λg(cni)pni(1(rni)vni)λf(cni)pn+1i=12(Δx)2kRi([Dp(Mn+1k+Mni)DhNn+1k(vn+1kvn+1i)]pn+1k)12(Δx)2kRi(Dp(Mn+1k+Mni)pn+1iDhNni(vn+1kvn+1i)pni)+λg(cni)pni(1(rni)vni)λf(cni)pn+1i

    where spn+1i+1=pn+1i+1pn+1i+1+qn+1i+1, spn+1i1=pn+1i1pn+1i1+qn+1i1, spni=pnipni+qni,

    Mn+1i+1=κn+1i+11+spn+1i+1vn+1i+1, Mn+1i1=κn+1i11+spn+1i1vn+1i1, Mni=κni1+spnivni,

    Nn+1i+1=κn+1i+1vn+1i+11+vn+1i+1pn+1i+1, Nn+1i1=κn+1i1vn+1i11+vn+1i1pn+1i1 , Nni=κnivni1+vnipni.

    The index set Ri={i1,i+1} points at the direct neighbors of the space node xi. The discretization of the equation describing the proliferating cancer cells can thus be written as:

    Appn+1=pn+ϑnp (4.2)

    The tridiagonal matrix Ap has size (k+1)×(k+1), pn+1 is the vector containing the values of p for the space node k+1 at the (n+1)th time level. The vector ϑnp contains the entries Dtλg(cni)pni(1rnivn+1i) for i=1,2,3,...,k. Similarly, the equation for the quiescent cells read:

    qn+1iqniDt=Dq2(Δx)2kRi(Φq(χn+1k)+Φq(χni))(qn+1kqn+1i)+f(cni)pnih(cni)qn+1i (4.3)

    where χn+1i+1=κn+1i+11+sqn+1i+1vn+1i+1, χn+1i1=κn+1i11+sqn+1i1vn+1i1, χni=κni1+sqnivni.

    Equation (4.3) can also be written as

    Aqqn+1=qni+ϑnq (4.4)

    The tridiagonal matrix Aq has size (k+1)×(k+1), qn+1 is the vector containing the values of q for the space node k+1 at the (n+1)th time level. The vector ϑnq contains the entries Dtf(cni)pni for i=1,2,3,...,k. From the necrotic cell density and the nutrient concentration we obtain:

    un+1i=uni+Dtδuqni.
    cn+1icniDt=Dc(Δx)2ki=1[(cn+1i+1cn+1i)+(cn+1i1cn+1i)]+k1(1α(pni+qni))k1cn+1ik2pnicn+1i
    Accn+1i=cni+ϑnc (4.5)

    where ϑnc is the vector containing Δtk2c0(1α(pni+qni)) for i=1,2,3,...,k.

    The following initial conditions were used in the numerical simulations for proliferating cell, quiescent cell, ECM, necrotic cell, a proteolytic enzyme, and contractivity densities respectively.

    p0(x)=exp(x20.01),x[0,1]q0(x)=0.3p0(x),x[0,1]v0(x)=1exp(x20.01),x[0,1]u0(x)=0,x[0,1]y0(x)=15fγ(80x,3,8)+0.2,x[0,1]κ0(x)=2y0(x),x[0,1].

    As mentioned above, early detection of tumor cells is difficult therefore, we initially assume a small lump of tumor cells have penetrated the tissue.

    Table 1 shows the parameter description, range, and values used in the model.

    Table 1.  Table of parameters used in the model.
    Parameter Description Value Range Source
    Dc (Diffsion coefficient of nutrient) 1 1022 [15]
    Dp (Diffusion coefficient of proliferating cells) 102 104102 [14]
    Dq (Diffusion coefficient of quiescent cells) 102 104102 [14]
    Dh (Haptotaxis coefficient) 1 10310 [6,14]
    δv (Rate of ECM degradation) 15 120 [6]
    a (decay rate) 3 15 [6]
    μv (Rate of ECM proliferation) 0.3 0.152.5 [3,6]
    c0 (Initial concentration of nutrients) 1 01 [15]

     | Show Table
    DownLoad: CSV

    In this section, we show the numerical solutions for the densities of three cell layers (proliferating, quiescent, necrotic), ECM, and contractivity function in space and specific intervals of time (t=1,2,5,10,...50). Figure 1 (A-L), Figure 2 (A-F), and Figure 3 (A-C) show qualitative solutions without delay (κ=1, solid lines) and with a delay (τ=5 dashed lines) for: η=0.01,λ=2; η=0.01,λ=0.02: and η=0.005,λ=1.5 respectively.

    Figure 1.  Evolution of proliferating tumor cell density (blue), quiescent cell density (pink), necrotic cell density (red), ECM fiber density (green), and contractivity function (black) in the pure macroscopic setting (solid lines) and with the delay of τ=5 (dashed line), η=0.01, M=1, λ=2 for equations (3.1) and (2.4).
    Figure 2.  Evolution of proliferating tumor cell density (blue), quiescent cell density (pink), necrotic cell density (red), ECM fiber density (green), and contractivity function (black) in the pure macroscopic setting (solid lines) and with the delay of τ=5 (dashed line), η=0.01, M=1, λ=0.02 for equations (3.1) and (2.4).
    Figure 3.  Evolution of proliferating tumor cell density (blue), quiescent cell density (pink), necrotic cell density (red), ECM fiber density (green), and contractivity function (black) in the pure macroscopic setting (solid lines) and with the delay of τ=5 (dashed line), η=0.005, M=1, λ=1.5 for equations (3.1) and (2.4).

    As earlier mentioned, biochemical reaction at the subcellular level takes place on a short time scale whiles its effects on the motile behavior of tumor cells take a much longer time to happen. We coupled the two-time scales by introducing the delay term in the contractivity function which we modeled as an input to the nonlinear density-dependent diffusion function of proliferating and quiescent cells, as well as the haptotaxis coefficient of proliferating cells. As time progress, proliferating cells invade faster without the delay (κ=1) than with delay as expected. Therefore it can be concluded that subcellular events slow down the invasion of proliferating cells into the tissue. Furthermore, there are higher peaks at the leading edge for proliferating cells (this may be due to the cells encountering undegraded tissue or bone).

    The density of necrotic cells builds up gradually at the center (x=0) as time increases. Accounting for the subcellular events shows that they seem to have no effect on the evolution of quiescent and necrotic cell densities. The invasion of ECM fiber is as expected. From Figures 1, 2 and 3, we notice that the solution has not reached the final time. Therefore, there is negligible influence of the boundary conditions on the invasion patterns.

    The study makes the following contributions to existing mathematical models on tumor invasion of host tissues: Our numerical solution for the system of equations showed the three distinct layers of proliferating, quiescent and necrotic cell densities. The subcellular, microscopic scale and macroscale events have varying spatial and time scales, we coupled the two-time scales using a contractivity function. As a result, we were able to determine the effects of biochemical dynamics at the subcellular level on the strength and motility of tumor cells. The model accounts for contact inhibition of proliferating and quiescent cells with the tumor micro-environments.

    The authors declare there is no conflicts of interest in this paper.



    [1] H. Jan, R. Kathy, W. Wendy, et al. A guide to cancer drug development and regulation, Astra Zeneca, 2006.
    [2] F. Cornelis, O. Saut, P. Comsille, et al. In vivo mathematical modeling of tumor growth from imaging data: Soon to come in the future? Diagn. Interv. Imag., 94 (2013), 593-600.
    [3] G. Meral, C. Stinner and C. Surulescu, On a multiscale model involving cell contractivity and its effects on tumor invasion, Discrete and continuous dynamical systems series B, 20 (2014), 189-213. doi: 10.3934/dcdsb.2015.20.189
    [4] J. Kelkel and C. Surulescu, On some models for cancer cell migration through tissue networks, Math. Biosci. Eng., 8 (2011), 575-589. doi: 10.3934/mbe.2011.8.575
    [5] B. Yue, Biology of the Extracellular Matrix: An Overview, J. Glaucoma, 23 (2014), 20-23. doi: 10.1097/IJG.0000000000000108
    [6] M. A. J. Chaplain, G. Lolas, Mathematical modelling of cancer invasion of tissue: Dynamic heterogeneity, Netw. Heterog. Media, 1 (2006), 399-439. doi: 10.3934/nhm.2006.1.399
    [7] N. Kolbem, A tumor invasion model for heterogeneous cancer cell populations: mathematical analysis and numerical methods, Ph.D thesis, der Johannes Gutenberg-Universität Mainz, 2017. Available from: https://https://publications.ub.uni-mainz.de/theses/volltexte/2018/100001842/pdf/100001842.pdf.
    [8] N. Sfakianakis, N. Kolbe, N. Hellmann, et al. 2016. A multiscale approach to the migration of cancer stem cells: mathematical modelling and simulations, B. Math. Biol., 79 (2017), 209-235. doi: 10.1007/s11538-016-0233-6
    [9] P. Domschke, D. Trucu, A. Gerisch, et al. Mathematical modelling of cancer invasion: Implications of cell adhesion variability for tumor infinitive growth patterns, J. Theor. Biol., 361 (2014), 41-60. doi: 10.1016/j.jtbi.2014.07.010
    [10] V. Andasari, A. Gerisch, G. Lolas, et al. Mathematical modelling of cancer cell invasion of tissue: biological insight from mathematical analysis and computational simulation, J. Math. Biol., 63 (2011), 141-171. doi: 10.1007/s00285-010-0369-1
    [11] P. Lu, T. Dumitru, L. Ping, et al. Multiscale Mathematical Model of Tumour Invasive Growth, B. Math. Biol., 79 (2017), 389-429. doi: 10.1007/s11538-016-0237-2
    [12] T. Dumitru, Multiscale Modelling of Cancer Invasion, 6th European Conference on Computational Mechanics (ECCM 6), 2018.
    [13] R. M. Sutherland, Cell and environment interaction in tumour microregions: the multicell spheroid model, Science, 240 (1988), 177-184. doi: 10.1126/science.2451290
    [14] N. Kolbe, M. Lukacova-Medvid ova, N. Sfakianakis, et al. Numerical simulation of a contractivity based multiscale cancer invasion model, Institute of Math., Johannes Gutenberg-Uni., Mainz, Germany, 2016.
    [15] J. A. Sherratt, M. J. A. Chaplain, A new mathematical model for avascular tumour growth, J. Math. Biol., 43 (2001), 291-312. doi: 10.1007/s002850100088
    [16] Y. Tao, C. Cui, A density-dependent chemotaxis-haptotaxis system modeling cancer invasion, J. Math. Anal. Appl., 367 (2010), 612-624. doi: 10.1016/j.jmaa.2010.02.015
    [17] R. A. Gatenby, E. T. Gawlinski, A reaction-diffusion model of cancer invasion, Cancer Res., 56 (1996), 5745-5753.
    [18] P. Domschke, D. Trucu, A. Gerisch, et al. Mathematical modelling of cancer invasion: Implications of cell adhesion variability for tumour infiltrative growth patterns, J. Theor. Biol., 361 (2014), 40-61. doi: 10.1016/j.jtbi.2014.07.010
    [19] A. Zhigun, C. Surulescu, A. Hunt, Global existence for a degenerate haptotaxis model of tumor invasion under the go-or-grow dichotomy hypothesis, preprint.
    [20] W. Mueller-Klieser, Method for determination of oxygen consumption rates and diffusion coefficients in multicellular spheroids, Biophys. J., 46 (1984), 343-348. doi: 10.1016/S0006-3495(84)84030-8
    [21] W. Mueller-Klieser, J. P. Freyer, R. M. Sutherland, Influence of glucose and oxygen supply conditions on the oxygenation of multicellular spheroids, Brit. J. Cancer, 53 (1986), 345-353. doi: 10.1038/bjc.1986.58
    [22] T. Roose, S. J. Chapman, P. K. Maini, Mathematical Models of Avascular Tumor Growth, Siam rev., 49 (2007), 179-208. doi: 10.1137/S0036144504446291
    [23] T. S. Deisboeck, Z. Wang, P. Macklin, et al. Multiscale Cancer Modeling, Annual review of biomedical engineering, 2010.
    [24] D. Hanahan and R. A. Weinberg, The Hallmarks of Cancer, Cell, 100 (2000), 57-70. doi: 10.1016/S0092-8674(00)81683-9
    [25] D. Hanahan, R. A. Weinberg, Hallmarks of cancer: the next generation, Cell, 144 (2011), 646-674. doi: 10.1016/j.cell.2011.02.013
    [26] B. Qian, J. W. Pollard, Macrophage Diversity Enhances Tumor Progression and Metastasis, Cell, 141 (2010), 39-51. doi: 10.1016/j.cell.2010.03.014
    [27] J. A. Joyce, J. W. Pollard, Microenvironmental regulation of metastasis, Nat. Rev. Cancer, 9 (2009), 239-252. doi: 10.1038/nrc2618
    [28] M. A. J. Chaplain and G. Lolas, Mathematical modelling of cancer cell invasion of tissue: The role of the urokinase plasminogen activation system, Math. Mod. Meth. Appl. S., 15 (2005), 1685-1734. doi: 10.1142/S0218202505000947
    [29] V. Andasari, R. T. Roper, M. H. Swat, et al. Integrating Intracellular Dynamics Using CompuCell3D and Bionetsolver: Applications to Multiscale Modelling of Cancer Cell Growth and Invasion, PLoS ONE, 7 (2012), e33726. doi: 10.1371/journal.pone.0033726
    [30] M. A. J. Chaplain and N. Deakin, Mathematical modeling of cancer invasion: the role of membrane-bound matrix metalloproteinases, Frontiers in oncology, 2013.
    [31] C. Walker, E. Mojares, A. D. R. Hernández, Role of Extracellular Matrix in Development and Cancer Progression, Int. J. Mol. Sci., 19 (2018), 3028.
    [32] Y. T. N. Edalgo, A. N. F. Versypt, Mathematical Modeling of Metastatic Cancer Migration through a Remodeling Extracellular Matrix, Processes, 6 (2018), 58.
    [33] P. J. Murray, C. M. Edwards, M. J. Tindall, et al. From a discrete to a continuum model of cell dynamics in one dimension, Phys. Rev. E, 80 (2009), 31912. doi: 10.1103/PhysRevE.80.031912
    [34] A. Gerischa, M. A. J. Chaplain, Robust numerical methods for taxis-diffusion-reaction systems: Applications to biomedical problems, Math. Comput. Model., 43 (2006), 49-75. doi: 10.1016/j.mcm.2004.05.016
  • This article has been cited by:

    1. Ashlee N. Ford Versypt, Multiscale modeling in disease, 2021, 27, 24523100, 100340, 10.1016/j.coisb.2021.05.001
    2. D. Andreucci, A.M. Bersani, E. Bersani, F.J. León Trujillo, S. Marconi, A 3D mathematical model of coupled stem cell-nutrient dynamics in myocardial regeneration therapy, 2022, 537, 00225193, 111023, 10.1016/j.jtbi.2022.111023
    3. Mahmoud B. A. Mansour, Hussien S. Hussien, Asmaa H. Abobakr, Numerical simulations of wave propagation in a stochastic partial differential equation model for tumor–immune interactions, 2022, 0, 1565-1339, 10.1515/ijnsns-2022-0026
    4. Zijing Ye, Shihe Xu, Xuemei Wei, Analysis of a free boundary problem for vascularized tumor growth with a time delay in the process of tumor regulating apoptosis, 2022, 7, 2473-6988, 19440, 10.3934/math.20221067
    5. Piyush Pratap Singh, Binoy Krishna Roy, Chaos and multistability behaviors in 4D dissipative cancer growth/decay model with unstable line of equilibria, 2022, 161, 09600779, 112312, 10.1016/j.chaos.2022.112312
    6. Nabeela Anwar, Iftikhar Ahmad, Adiqa Kausar Kiani, Muhammad Shoaib, Muhammad Asif Zahoor Raja, Intelligent solution predictive networks for non-linear tumor-immune delayed model, 2024, 27, 1025-5842, 1091, 10.1080/10255842.2023.2227751
    7. Ali Sadiq Alabdrabalnabi, Ishtiaq Ali, Stability analysis and simulations of tumor growth model based on system of reaction-diffusion equation in two-dimensions, 2024, 9, 2473-6988, 11560, 10.3934/math.2024567
    8. Chien-Han Yen, Yi-Chieh Lai, Kuo-An Wu, Morphological instability of solid tumors in a nutrient-deficient environment, 2023, 107, 2470-0045, 10.1103/PhysRevE.107.054405
    9. Bengisen Pekmen, Ummuhan Yirmili, Numerical and statistical approach on chemotaxis-haptotaxis model for cancer cell invasion of tissue, 2024, 4, 2767-8946, 195, 10.3934/mmc.2024017
    10. M. B. A. Mansour, Hussien S. Hussien, Asmaa H. Abobakr, Numerical simulations of wave patterns in a stochastic partial differential equation model of tumor growth and invasion, 2025, 0971-3514, 10.1007/s12591-025-00720-6
    11. Pablo Facundo Ríos Flores, ¿Afectos populistas para una democracia radical? Ernst Cassirer y los riesgos de la convocatoria populista a la potencia afectiva del mito, 2025, 81, 2386-5822, 219, 10.14422/pen.v81.i313.y2025.012
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(7296) PDF downloads(668) Cited by(11)

Figures and Tables

Figures(3)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog