A discrete in continuous mathematical model of cardiac progenitor cells formation and growth as spheroid clusters (Cardiospheres)
Ezio Di Costanzo, Alessandro Giacomello, Elisa Messina, Roberto Natalini, Giuseppe Pontrelli, Fabrizio Rossi, Robert Smits, Monika Twarogowska
AA discrete in continuous mathematical model ofcardiac progenitor cells formation and growth asspheroid clusters (Cardiospheres)
Ezio Di Costanzo ∗ , Alessandro Giacomello ∗∗ , Elisa Messina ∗∗ , Roberto Natalini ∗ , Giuseppe Pontrelli ∗ , Fabrizio Rossi ∗∗ , Robert Smits ∗∗∗ , Monika Twarogowska ∗∗ Istituto per le Applicazioni del Calcolo “M. Picone”Consiglio Nazionale delle RicercheVia dei Taurini 19 – 00185 Rome, Italy. ∗∗ Department of Molecular MedicinePasteur Institute Cenci–Bolognetti FoundationSapienza University of Rome, Viale del Policlinico – 00161 Rome, Italy. ∗∗∗
Department of Mathematical SciencesNew Mexico State UniversityLas Cruces, New Mexico, USA 88003.
Abstract
We propose a discrete in continuous mathematical model describing the in vitrogrowth process of biophsy-derived mammalian cardiac progenitor cells growing as clustersin the form of spheres (
Cardiospheres ). The approach is hybrid: discrete at cellularscale and continuous at molecular level. In the present model cells are subject to theself-organizing collective dynamics mechanism and, additionally, they can proliferateand differentiate, also depending on stochastic processes . The two latter processes aretriggered and regulated by chemical signals present in the environment. Numericalsimulations show the structure and the development of the clustered progenitors andare in a good agreement with the results obtained from in vitro experiments.
Keywords.
Mathematical Biology, differential equations, hybrid models, Poisson stochasticprocess, collective dynamics, cell movements, cellular signaling, chemotaxis, stem cells,Cardiospheres.
Mathematics Subject Classication.
E-mail: [email protected] , [email protected] , [email protected] , [email protected] , [email protected] , [email protected] , [email protected] , [email protected] a r X i v : . [ q - b i o . CB ] D ec Introduction
First described in 2004, the cardiac biopsy-derived progenitor cells, growing in vitro as“niche”-like microtissue, known as
Cardiospheres (CSps), represent a widely adopted platformtechnology holding great promise to realize a powerful cell therapy system, in vitro drugscreening and disease modeling opening new opportunities for myocardial repair (Messinaet al, 2004; Ruckdeschel Smith et al, 2007; Li et al, 2010). To achieve this goal, an optimalmanagement of qualitative and quantitative CSps growth-modifying factors in terms of cellproliferation, differentiation and interaction with the external environment has to be known.It requires a simulator of the features and the fate of the Cardiospheres-niche-system whichis based on new as well as already existing experimental biological data and on reliablemathematical models.Our aim is to model the formation and differentiation of the original 3D cardiac biopsy-derived heart progenitors in vitro system, that is the CSp. These structures are cellularspheroids with a central nucleus of less differentiated elements surrounded by outer layersof cells more committed toward different levels of cardio-vascular differentiation with otherproliferation features (Figure 1(a)). Their individual formation in vitro starts from a few(low differentiated and less specialized) progenitor cells which are cultured in a carefullychosen environment and under specific conditions (Messina et al, 2004; Chimenti et al, 2012),and have the potential to develop into a more cardiac-specific differentiation which, unlessderived from pre-natal hearts, are no longer able to achieve a terminal level of contractilecardiomyocytes. The major types of molecular processes that control cellular differentiationand proliferation involve various physical factors, nutrients and cell signaling. During themultiple stages of development, the cell size and shape changes dramatically, as a responseto signaling molecules (Brown et al, 2003; Agathocleous and Harris, 2013). This biologicalenvironment is similar to that of tumor spheroids as described in Wallace and Guo (2013).Neoplastic cells hold several growth strategies, such as genetic mutation of cell cycle inhibitors,neoangiogenesis, etc., and cannot be directly compared with those more limited of the adultcardiac progenitors. Nevertheless both of the 3D cellular growth mechanisms share somecommon features. CSps and neoplastic spheroids hold more than two functional phenotypes(proliferating, quiescent, differentiating, dying). The growth in a form of a 3D sphere limitsthe diffusion of the nutrients and oxygen and leads to the formation of a central core of deathcells as a consequence of necrosis, apoptosis and anoichisis (Figure 2). It is well known thatin tumors the fastest proliferating cells are present in the intermediate and external layersforming the so-called growing front. On the contrary, in small and medium CSps the highestrate of the proliferation is observed in the central area due to the hypoxic conditions, which arecloser to those corresponding to the pre-natal developing heart. However, while the spheresbecome larger the necrotic core appears as in tumors. In this paper we introduce a discrete incontinuous non-deterministic model to simulate the growth of a CSp aggregate. We treat thevarious biological constituents using different levels of description: on the cellular scale weadopt a discrete description, while on the molecular scale chemical signals are considered ascontinuous variables (Anderson et al, 2007; Preziosi and Tosin, 2009). This choice is justifiedby the relatively small number of cells involved in the formation of CSps, for which thehypothesis of the continuum does not hold, and by the much smaller size of the moleculeswith respect to them. Similar approaches have been adopted in analogous cell systems, forinstance in relation to the modeling of tumor growth (Enderling et al, 2012; Kim, 2013;2 a) (b)
Figure 1: (a) A CSp by immunofluorescence after 3 days: a cluster of immature stemcells (green) is localized in the central core of Cardiospheres, surrounded by differentiatedsupporting cells (red). Blue spots indicate cell nuclei, and cells inside the black area (nonlabeled) are to be intended in an intermediate level of differentiation, see also Section 2(source from Li et al, 2010, by courtesy of John Wiley & Sons). (b) An example of anumerical simulation of our mathematical model showing a CSp at 3 days. Green, blue andred colors mark respectively three increasing degrees of maturation of a single cell (see Section4 for further details).Figure 2: Histological preparation of CSps ultra-thin slices. Hemathoxilin-eosin stainingshows a nuclei-free central necrotic area more evident in larger CSps.3sborne et al, 2015; Chaplain, 2011; Fletcher et al, 2015, and references therein) or, recently,to describe the morphogenesis of the lateral line in the zebrafish embryo (Di Costanzo et al,2015), and to study the behavior of the collective motion of cells under chemotaxis andalignment effect (Di Costanzo and Natalini, 2015). With particular regard to the stem cellcultures, in Wu et al (2014) a continuous diffusion–reaction mathematical model is proposedto investigate oxygen transport and distribution in embryonic stem cell aggregates, also incomparison with experimental data.Our proposed mathematical model relies on a hybrid description: each cell is describedseparately by a system of ordinary differential equations, which take into account the self–organizing interactions between cells, through typical terms of collective dynamics andchemotaxis effects (see for instance the reviews by Vicsek and Zafeiris, 2012; M´ehes andVicsek, 2014); at a macroscopic level, by reaction-diffusion equations, we describe theconcentrations of oxygen, as a representative of nutrients essential for cells metabolism, andof a growth factor TGF- β (Transforming-Growth-factor- β ), as a chemoattractant and mainregulator of the cell differentiation, secreted by the cells and also present in the environment.Additionally, to model the formation of the Cardiosphere, we introduce two processes, the cellproliferation and differentiation, which lead to the growth and maturation of an initial cellcluster. Moreover, we distinguish three differentiation levels of cell maturation, which alsoinfluence the two aforesaid processes. Moreover they are regulated by the local concentrationsof oxygen and TGF- β . Another crucial feature of our model is the stochasticity. It ismotivated by the so called noise in biological systems and the averaging of cells characteristics.The former is mainly due to the genetic diversity among individual cells, random collisionsand thermal fluctuations in chemical reactions and randomness of various external factors,see Tsimring (2014). The latter concerns the simplifications introduced in our model, suchas the uniform size of cells, equal maturation time or signal sensing. The randomness of thebiological phenomena in our model is present in the proliferation and differentiation processes.When the necessary conditions are satisfied cells divide or pass to a higher differentiated levelwith some probability, which increases if the conditions are favorable.Using the proposed mathematical model we have simulated the formation of aCardiosphere in a two–dimensional setting. After a careful sensitivity analysis of the modelparameters we were able to obtain numerical results comparable with the existing biologicalexperiment: see again Figure 1, in which the structure and composition of the two images,real experiment (a) and numerical simulation (b), are in good agreement. In order to meet theneed of biological research and possible therapeutic applications (including the achievement ofa more advanced differentiation level) we have analyzed the development of the Cardiosphereat two typical experimental concentration of oxygen.The paper is organized as follows: in Section 2 we present the general setting of theproblem and the mathematical framework used throughout the paper. In Section 3 wedescribe in detail the mathematical model for the cells and nutrients dynamics, togetherwith the stochastic mechanisms of cell proliferation and differentiation. Section 4 containsthe numerical results and relative biological interpretations. Finally, Section 5 is devoted tothe conclusions and possible future perspectives.4 Formulation of the problem
The in vitro 3D CSp formation starts from a low-density seeding of the explants-biopsy-derived cells in a defined culture conditions. First, adhered to the culture dish bottomthe individual stem cells start to aggregate in small clusters of very few elements. Reachingcertain dimensions they detach and grow immersed in a liquid environment and in a controlledgas atmosphere, which contain substances essential for their metabolism. In our model wehave reduced to two the biological and physical variables affecting the CSps growth anddifferentiation: the oxygen, representing nutrients and gases, and the TGF- β which is presentin the culture medium and is produced by the cells themselves. The latter represents a keymodulator of the main pathway involved in the induction of the non-terminal differentiationlevel achieved in our basic experimental conditions (Forte et al, 2012). We assume that duringthe formation of a CSp, the cells are subject to two basic cellular processes, the proliferationand the differentiation, which depend, with some probability, on the type of a cell and on thelocal availability of oxygen. Moreover, we consider three different increasing differentiationlevels for living cells and, separately, dead cells. The complex differentiation mechanism istriggered by a large enough concentration of TGF- β , is hindered by the presence in the closestproximity of a large number of cells of the same type and is inhibited when the concentrationof oxygen drops below a threshold value (see Section 3.2).To describe the formation of CSps in the next section we propose a hybrid model, inwhich cells are considered as discrete entities while oxygen and TGF- β concentrations arecontinuum variables. We assume that in the absence of specific signals the cell dynamics isdue to attractive-repulsive forces and friction effect, taking place isotropically. The TGF- β ,produced by the cells and also present in the external medium, is a chemoattractant forcell which moves towards its higher concentrations. For computational simplicity a two-dimensional description is adopted. In a Cartesian coordinate system, let us consider a rectangular domain Ω, which in our modelrepresents a limited portion of the culture (containing many CSps) with a single cluster ofstem cells in its center (Figure 3). For this reason in the following we can assume periodicboundary condition in relation to the equations of the model. At the initial time, N cellsat the same maturation level (Section 3.2) and radius R (for simplicity all cells are assumedto have the same constant radius) are placed randomly at positions X i , i = 1 , ,...,N at thecenter of the domain. In the following, due to proliferation, the number of cells at the time t will be denoted by N ( t ). In the model we denote by c ( x ,t ) and S ( x ,t ) the concentrationsof oxygen and TGF- β respectively, and in the numerical simulations these variables arediscretized on a uniform spatial mesh. In general, the cell position X i does not coincide witha mesh grid point, and we have to “approximate” the value of c ( X i ,t ). In general, a cell”feels” the presence of chemical substance (or a signal) h ( X i ,t ) not only at the center X i ,but also in its proximity. To model this graded influence, a weighted average operator isdefined: F ( h ( X i ,t )) = 1 W (cid:90) B ( X i , ¯ R ) w i ( x ) h ( x ,t ) d x , (1)5igure 3: Experimental picture of a culture of CSps. The dotted red line indicates the portionchosen as possible domain Ω for a single CSp in our mathematical model.where: W := (cid:90) B ( X i , ¯ R ) w i ( x ) d x , (2)and a possible choice for the weight function can be a truncated Gaussian w i ( x ) := (cid:18) −|| x − X i || log 2¯ R (cid:19) − , if || x − X i || ≤ ¯ R, , otherwise , which models a sensing intensity decreasing on the distance from the cell center. Here ¯ R represents an influence radius, possibly greater than the cell radius R (Table 2), and B ( X i ,R ) := { x : || x − X i || ≤ R } is the ball centered in X i with radius R . In the following the characteristic function over theball B ( X i ,R ) is defined as: χ B ( X i,R ) := (cid:26) , if x ∈ B ( X i ,R ) , , otherwise . In this section we present in detail the discrete in continuous model for the formation of aCardiosphere. Firstly, we present the equations describing the movement of cells and theforces acting on them. Then we focus on the proliferation and differentiation processes. Thebasic mechanisms of the phenomena are explained and the reaction-diffusion systems for thetime evolution of the oxygen and TGF- β concentrations are given. Motion of cells is the result of mechanical forces acting on them and of their reciprocalinteractions. Due to the finite volume of cells there is an internal pressure between them6hich leads to repulsion while being compressed. On the other hand, there are attractivebonds between cells which keep them in contact. We also assume that cells are driven by achemotactic signal, that is they move toward the higher gradient of the TGF- β concentration,and are slowed down because of the viscosity of the medium in which they are immersed(Rubinstein et al, 2009; Fournier et al, 2010; Bayly et al, 2012).The position of the i -th cell X i follows Newton’s law and is described by the second orderordinary differential equation¨ X i = (cid:88) j : X j ∈ B ( X i ,R ) \{ X i } K ( r ij ) − µ ˙ X i + α F ( ∇ S ( X i ,t )) , (3)where dot denotes time derivative, µ is a friction coefficient for unit mass and K ( r ij ) is theattraction-repulsion force for unit mass between i-th and j-th cell, defined as K ( r ij ) := − k (cid:18) || r ij || − R (cid:19) r ij || r ij || , || r ij || ≤ R ,k ( || r ij || − R ) r ij || r ij || , R < || r ij || ≤ R . (4)The parameters k , k are positive constants, r ij = X j − X i is the distance between cells, R = 2 R and R (to be chosen later, see Table 2) are the radii of repulsion and attractionrespectively. Note that equation (4) gives a repulsion in the form 1 /r , while (4) is a linearattractive elastic force. Similar terms can be found in the model proposed in D’Orsogna et al(2006); Joie et al (2013). Equation (3) is augmented with the initial conditions X i (0) = X i , ˙X i (0) = V i = . Some of the present authors has been investigating the proliferation and differentiationdynamics of cells forming CSps from the original description of the method of isolation andexpansion of human cardiac progenitors from heart biopsy (Messina et al, 2004). Confocalanalysis of the BrdU labeled cells which incorporate the base analogue within the replicatingDNA, showed the positive fluorescence labeled cells located in the inner of the growingspheres. No or a weak signal was present in the external layers of more differentiated cells(that cannot replicate at the same speed). This kind of layered growth implies that theproliferation in the central area depends also on the CSps size and the core-cells may becomequiescent or die by necrosis or apoptosis.It is worth noting that tumor spheroid growth depends (at least partially) on the samefactors, however, in this case they show a more pronounced variety of final volumes for anygiven cell line. We can distinguish an early phase of global undifferentiated mass, a secondstage of a shell of proliferating cells around a central area of quiescent elements and thelast stage in which an inner necrotic/apoptotic core is enveloped by the previous two. Theexternal layer always represents the growing front (Kim, 2013).Stem cells are characterized by different levels of maturation (differentiation). To modelthis feature we mark the i -cell with a label ϕ i ( t ) = 1 , , ,d , called its “state”. The state ϕ = d indicates dead cells or the space occupied by necrotic mass (due to adverse biological7onditions, as oxygen deficiency for a sufficiently long period of time). Note that apoptosis,which is a programmed cells death, occurs over a longer timescale and is not considered here.In particular, we introduce a threshold concentration c d such that whenever F ( c ( X i ,t )) > c d cell survives, otherwise it dies. Such cell is marked as dead ( ϕ = d ) and does not take part infurther evolution of the CSp, that is it neither proliferates nor differentiates. However, deadcells do occupy a volume and exert attractive-repulsive forces on other cells. A stochasticprocess introduced in the process of proliferation and differentiation models also the presenceof quiescent cells. Such cells neither proliferate nor differentiate but, in contrast to deadcells, this state results from environment conditions and is only temporary. At a later timea quiescent cell can restart its functions.For a living cell, ϕ = 1 denotes the lowest degree of maturation (the least differentiatedcells), while ϕ = 2 , cell cycle time T c , i.e. a cell maturationtime after which the cell may proliferate and/or differentiate with some probability. In the model we assume that, with some probability and when environment conditions aresatisfied, from time t to time t + ∆ t a cell can change its state to a higher one by a unityaccording to Table 1, or can die (as described in Section 3.2.2). ϕ i ( t + ∆ t ) = 1 ϕ i ( t + ∆ t ) = 2 ϕ ( t + ∆ t ) = 3 ϕ i ( t ) = 1 (cid:88) (cid:88) ϕ i ( t ) = 2 (cid:88) (cid:88) ϕ i ( t ) = 3 (cid:88) Table 1: Possible differentiation levels reached by a cell. The check marks the admissibletransition at time t + ∆ t for each cell state on the first column.More precisely, the differentiation is triggered by chemical signals, such as TGF- β presentin the extracellular environment and secreted by cells themselves. Moreover, it can beinhibited due to the presence, in the surrounding area, of cells of the same type, as wecan found in similar cell systems (Itoh and Chitnis, 2001; Matsuda and Chitnis, 2010) andmathematical models (e.g. Di Costanzo et al, 2015). We denote by Γ i an inhibitor indicatorfunction, which depends on the number n i of cells in the neighborhood of the i-th cell and isdefined asΓ i ( t ) = (cid:26) , if n i ≤ ¯ n, , if n i > ¯ n, n i := card { j : X j ∈ B ( X i ,R ) } , i = 1 ,...,N. For the threshold ¯ n we set in the following the value ¯ n = 4 for a two dimensional case.We define a differentiation threshold q i ( t ) := σ ( ϕ i ) · F ( S ( X i ,t )) · Γ i ( t ) S max , i = 1 ,...,N, (5)for the cell i , with σ ( ϕ i ) ( ϕ i = 1 ,
2) positive constants, and S max the experimental maximumvalue of TGF- β (Table 2). In addition, we assume that differentiation occurs under the8ondition that oxygen concentration level is above a given threshold value ¯ c ( ϕ i ), alsodepending on the cell states ϕ i = 1 , ϕ i ( t ) can be modeled as a non-homogeneous Poisson processwith variable intensity q i ( t ) (see for example Anderson et al, 2015 and references therein, aswell as Kusher and Dupuis, 1992). In particular, the probability of a differentiation in a timeinterval ∆ t is given by P ( ϕ i ( t + ∆ t ) − ϕ i ( t ) ≥ (cid:117) q i ( t )∆ t, (6)the expected number of differentiations in the interval ( t,t + ∆ t ) is E (number of differentiations in ( t,t + ∆ t )) = q i ( t )∆ t, (7)therefore in the time interval (0 ,t ) results E (number of differentiations in (0 ,t )) = (cid:90) t q i ( s ) ds. (8)Finally, we write the stochastic equation dϕ i = d P (cid:18)(cid:90) t q i ( s ) ds (cid:19) , i = 1 ,...,N, (9) P ( (cid:82) t q i ( s ) ds ) being a inhomogeneous Poisson process with intensity q i ( t ).In practice, if F ( c ( X i ,t )) ≥ ¯ c ( ϕ i ), starting from (9), we compute the threshold (5), and wecan adopt the following stochastic rule: cell types ϕ i = 1 and ϕ i = 2 switch irreversibly theirstate according to ϕ i ( t + ∆ t ) = ϕ i ( t ) + (cid:26) , if ¯ q i < q i ∆ t, , otherwise , i = 1 ,...,N, (10)where ¯ q i is a random number between 0 and 1 (see also Section 4). We introduce the proliferation thresholds p i ( t ) of the i-th cell at time t in the following form p i ( t ) = g ( F ( c )) , if ϕ i ( t ) = 1 ,g ( F ( c )) , if ϕ i ( t ) = 2 ,g ( F ( c )) , if ϕ i ( t ) = 3 . (11)Functions g ,g ,g describe the dependence of the rate of proliferation on the concentration ofoxygen. Least differentiated cells, i.e. those at state ϕ = 1, work under anaerobic conditions.Their proliferative ability increases at low levels of oxygen with an assumed gaussian-typedependence. When c is too high, cells change their metabolism and proliferate at a lowerrate. Cells at states ϕ = 2 , ϕ = 1. As a consequence, a possiblechoice for functions g ,g ,g , used in the numerical simulations, is shown in Figure 4. If, for9ertain interval of time, the concentration of oxygen is below a given threshold c d , then cellsdo not have enough resources to live and die.To model stochasticity in proliferation, we introduce a new state variable φ i ( t ) denotingthe number of (from i − th ) generated cells in (0 ,t ). Similarly to the differentiation (Section3.2.1), φ i ( t ) satisfies a non-homogeneous Poisson process: dφ i = d P (cid:18)(cid:90) t p i ( s ) ds (cid:19) , i = 1 ,...,N, (12)with variable intensity p i ( t ) given by equation (11). P r o li f e r a t i on t h r e s ho l d s Cell type 1Cell type 2Cell type 3
Figure 4: Proliferation thresholds as a function of oxygen concentration (equation (11)).If F ( c ( X i )) > c d , we compute the proliferation threshold of the i-th cell p i through equation(11). The following computational rule of proliferation is applied: φ i ( t + ∆ t ) = φ i ( t ) + (cid:26) , if ¯ p i < p i ∆ t, , otherwise , i = 1 ,...,N, (13)where ¯ p i is a random number between 0 and 1 (see also Section 4).When a cell proliferates, a newborn ( N + 1)-th cell appears and is appended to the vectorof coordinates X i , i = 1 , ,...,N . Its position X N +1 is generated randomly within the distance r from the mother cell and with same velocity. More precisely, in polar coordinates wechoose randomly two numbers: the angle θ ∈ [0 , π ] and the distance between the two centers r ∈ [ R/ ,R ]. 10 .3 Dynamics of oxygen and TGF- β signal The concentration of oxygen c = c ( x ,t ) is a continuous variable and its time evolution isgoverned by a reaction-diffusion equation: ∂ t c = ∇ · ( D c ∇ c ) − N (cid:88) i =1 λ ( ϕ i ) c γ +1 k ( ϕ i ) + c γ +1 χ B ( X i,R + H ¯ B ( c − c ) ,c ( x ,
0) = c , in Ω , periodic boundary conditions, on ∂ Ω , (14)with c the environmental (constant) concentration. The last term in equation (14) modelsthe two-dimensional representation of the real 3D experiment: the supply of oxygen comingfrom the environment (upper surface of the culture) is recovered as a source term, modeledas proportional to the difference c − c with rate H , and spatially modulated by a well-shapedfunction over the CSp geometry, in the form¯ B ( x ) := exp (cid:34) ζ (cid:18) r ( x ) R CSp (cid:19) (cid:35) − ζ ) − , if r ≤ R CSp , , otherwise . (15)The constant ζ is a constant geometrical parameter, r ( x ) is the distance from the center ofthe sphere, and R CSp is the mean radius of CSp at time t . The shaped function (15) considersthe graded exposition to oxygen of cells at different depth. In the equation (14) we assumethat oxygen diffuses in the environment and is consumed by cells with a rate depending onthe metabolism of the cell (i.e. on its state ϕ i ), and increases with the oxygen availability( ∝ c γ , γ >
0) according to the saturation Michaelis-Menten-like law (Wu et al, 2014). A typicalvalue γ = 1 / S ( x ,t ) is governed by the similar IBV problem ∂ t S = ∇ · ( D S ∇ S ) + N (cid:88) i =1 ξ ( ϕ i ) χ B ( X i,R − ηS,S ( x ,
0) = S , in Ω , periodic boundary conditions, on ∂ Ω , (16)where ξ is the (state dependent) S -release rate, η is the molecular degradation rate and D S ( x )is a diffusion coefficient.Compared with the production term, a possible S -consumption termis extremely low and has been neglected.In addition, the diffusivity coefficient D (say D c in equation (14), D S in (16)) varies inrelation with the change of cell density, which we measured in a reference volume. Due tothe presence of a large number of cells, D is significantly reduced at the center of the CSp.Diffusion can be considered as in a porous medium, and the effective diffusivity D ( x ) can be11efined as: D ( x ) := D max ρ ¯ A ( x ) , where D max is the constant unperturbed diffusion, ρ is a compactness parameter and ¯ A ( x )is the occupancy fraction in a proper control volume containing x . In this section we present numerical simulations of a typical formation and growth of CSps,under typical conditions (Messina et al, 2004; Chimenti et al, 2012), using the previouslydescribed mathematical model.For convenience we summarize here our hybrid system of equations: ¨ X i = (cid:88) j : X j ∈ B ( X i ,R ) \{ X i } K ( r ij ) − µ ˙ X i + α F ( ∇ S ( X i )) ,dϕ i = d P (cid:18)(cid:90) t q i ( s ) ds (cid:19) ,dφ i = d P (cid:18)(cid:90) t p i ( s ) ds (cid:19) ,∂ t c = ∇ · ( D c ∇ c ) − N (cid:88) i =1 λ ( ϕ i ) c γ +1 k ( ϕ i ) + c γ +1 χ B ( X i,R + H ¯ B ( c − c ) ,∂ t S = ∇ · ( D S ∇ S ) + N (cid:88) i =1 ξ ( ϕ i ) χ B ( X i,R − ηS, (17)where K ( r ij ) and F ( ∇ S ( X i )) are given respectively in (4) and (1), while q i and p i given by(5) and (11). Initial and boundary conditions are given by: X i (0) = X i , ˙ X i (0) = ,N (0) = N , ϕ i (0) = 1 , ∀ i = 1 ,...,N ,c ( x ,
0) = c , periodic boundary conditions on ∂ Ω ,S ( x ,
0) = S , periodic boundary conditions on ∂ Ω . Initially, we consider a cluster of N = 15 cells at the first stage of differentiation ( ϕ i = 1, ∀ i ) placed at the center of a square domain Ω = [0 , × [0 , µ m ). Such setting modelsthe initiation of the numerical simulation at the time that corresponds to about 24 hours ofthe real experiment.Despite the low geometrical dimension and a number of simplifying assumptions, thelarge set of biological parameters influences the problem and their complete characterizationremains a difficult task. Each parameter is interconnected with the others and strongly affectsthe simulation. Some of them are taken from previously published results or literature, others12re chosen in a compatible and consistent range, and the remaining ones are varied to analyzethe sensitivity of the system. The complete set of biological parameters is given in Table 2. All the equations are discretized by finite differences. The second order equation (17) isreduced to a system of first order equations and solved by explicit Euler in time. In equations(17) , the diffusion terms are discretized by a standard centered difference in space andintegrated implicitly in time. The nonlinear reaction terms are treated explicitly in time.Finally, in equation (17) we have eliminated the stiff term − ηS using a classical exponentialtransformation.The square domain size is chosen sufficiently large, 225 µ m per side, so that the growingCSp does not reach the boundary in the typical time of observation. We consider an uniformmesh with grid spacing ∆ x = ∆ y = 3 . µ m, which guarantees a sufficient discretization overeach cell (having a diameter of 15 µ m, see Table 2) and a reasonable accuracy.The time step ∆ t has been fixed as the maximum value to ensure stability and non-negativity of the scheme (max∆ t = 0 .
02 h). For the same reason, when the oxygenconcentration drops below a given threshold, it is necessary to reduce the time step to∆ t = 0 .
001 h.
Mathematical models of biological phenomena are subject to sources of uncertainty, includingerrors of measurement, natural intrinsic variability of the system, absence of information andpoor or partial understanding of the driving forces and mechanisms. This uncertainty imposesa limit on our confidence in the response or output of the model. One of the challengesis to characterize the model parameters that are not identified experimentally, and theseparameters are marked as calibrated in Table 2.In this section we present the calibration procedure in order to obtain as reliable resultsas possible. To achieve this goal, we studied the influence of the parameters on the modeldynamics through a local sensitivity analysis (Saltelli et al, 2008; Clarelli et al, 2015). Thisapproach determines a degree of dependency between input parameters (one at a time) andthe results of simulations. Although our analysis does not take into account the presenceof interactions between parameters, as in the global sensitivity analysis, it can give usefulinformation for further exploration of the parameter space.We measure the sensitivity of the model to the variation of a positive quantity Ψ withrespect to a reference parameter p by defining a sensitivity index SV ≥ SV := | ¯Ψ( p ± ε ) − ¯Ψ( p ) | / ¯Ψ( p ) ε/p , (18)where ¯Ψ is the average value of Ψ over a large number of runs (100 in our analysis) and itaccounts for the stochasticity presents in our model. Typically, ε is a small deviation over p ( ε = 0 . p , i.e. 5% variation). In Table 3 we report the sensitivity index computed withrespect to two observed variables: the total number of cells N , and the diameter of the CSp,both at the final simulated time of 72 hours. The smallness of the index give us confidenceof the robustness of the model in the response of the system in the presence of uncertainty.13 a b l e : E s t i m a t e s o f ph y s i c a l a ndb i o l og i c a l p a r a m e t e r v a l u e s . W h e n e v e r p o ss i b l e , t h e v a l u e s a r e t a k e n f r o m li t e r a t u r e o r e x p e r i m e n t s , t h e r e m a i n i n ga r ec a li b r a t e d t o b ec o n s i s t e n t w i t h o t h e r ( c . w . o . ) . T h e i r i nflu e n ce o n t h e s y s t e m i ss t ud i e d t h r o u g h a l o c a l s e n s i t i v i t y a n a l y s i s ( S ec t i o n . ) . W h e n a p a r a m e t e r a d m i t s a r a n g e o f p o ss i b l e v a l u e s , t h e u s e d o n e i n t h e nu m e r i c a l t e s t s i s pu t i nb r a c k e t s . P a r a m e t e r D e fin i t i o n E s t i m a t e d v a l u e o rr a n g e ( u s e d v a l u e s ) S o u r ce R ce ll r a d i u s . µ m W u e t a l ( ) ¯ R d e t ec t i o n r a d i u s o f c h e m i c a l s R b i o l og i c a l a ss u m p t i o n R r a d i u s o f a c t i o n o f r e pu l s i o nb e t w ee n ce ll s R b i o l og i c a l a ss u m p t i o n R r a d i u s o f a c t i o n o f a dh e s i o nb e t w ee n ce ll s . R b i o l og i c a l a ss u m p t i o n R r a d i u s o f p r o du c t i o n / d e g r a d a t i o n o f c h e m i c a l s R b i o l og i c a l a ss u m p t i o n R d e t ec t i o n r a d i u s f o r t h e d i ff e r e n t i a t i o n i nh i b i t i o n R b i o l og i c a l a ss u m p t i o n α c o e ffi c i e n t o f c h e m o t a c t i ce ff ec t p e r un i t m a ss µ m h − p g − c a li b r a t e d , c . w . o . k c o e ffi c i e n t o f r e pu l s i o np e r un i t m a ss µ m h − c a li b r a t e d , c . w . o . k e l a s t i cc o n s t a n t p e r un i t m a ss . × –1 . × ( . × ) h − B e ll e t a l ( ) µ f r i c t i o n c o e ffi c i e n t p e r un i t m a ss . × – ( . × ) h − R ub i n s t e i n e t a l ( ) D m a x c o xy g e nd i ff u s i o n c o e ffi c i e n t . . ( . ) × µ m h − W u e t a l ( ) D m a x S T G F - β d i ff u s i o n c o e ffi c i e n t . × µ m s − A I S P r o j ec t( ) ξ ( ϕ , , ) c o e ffi c i e n t o f p r o du c t i o n o f T G F - β . × − –1 . × − ( . × − ) p g µ m − h − C h i m e n t i e t a l ( ) S m a x m a x i m u m c o n ce n t r a t i o n o f T G F - β × − p g µ m − C h i m e n t i e t a l ( ) η d e g r a d a t i o n c o n s t a n t o f T G F - β . . ( . ) h − W a k e fi e l d e t a l ( ) σ ( ϕ = , ) d i ff e r e n t i a t i o n c o n s t a n t ( i f ϕ = ) , . ( i f ϕ = )( n o nd i m . ) c a li b r a t e d , c . w . o . λ ( ϕ = , , ) o xy g e n c o n s u m p t i o n c o n s t a n t . × − –2 . × − ( . × − i f ϕ = , . × − i f ϕ = , . × − i f ϕ = ) p g µ m − h − W u e t a l ( ) k ( ϕ = , , ) M i c h a e li s - M e n t e n o xy g e n c o n s t a n t . × − p g µ m − W u e t a l ( ) ¯ c ( ϕ = , ) o xy g e n c o n ce n t r a t i o n t h r e s h o l d v a l u e f o r d i ff . × − ( i f ϕ = ) , × − ( i f ϕ = ) p g µ m − c a li b r a t e d , c . w . o . c d o xy g e n m i n i m u m c o n ce n t r a t i o n f o r li f e ( n ec r o s i s ) . × − p g µ m − W u e t a l ( ) T a b l e : c o n t i n u a t i o n i n t h e n e x t pa ge a b l e : c o n t i n u a t i o n f r o m t h e p r ev i o u s pa ge P a r a m e t e r D e fin i t i o n E s t i m a t e d v a l u e o rr a n g e ( u s e d v a l u e s ) S o u r ce c e n v i r o n m e n t a l c o n ce n t r a t i o n f o r o xy g e n . × − p g µ m − ( i f O a t % ) , . × − p g µ m − ( i f O a t % ) W u e t a l ( ) S e n v i r o n m e n t a l c o n ce n t r a t i o n f o r S . × − p g µ m − T h e r m o S c i e n t i fi c ( ) T c ce ll c y c l e t i m e h ( i f O a t % ) , h ( i f O a t % ) M e ss i n a e t a l ( ) H s o u r ce r a t e h − c a li b r a t e d , c . w . o . ρ c o m p a c t n e ss p a r a m e t e r × ( n o nd i m . ) c a li b r a t e d , c . w . o . SV as defined in equation (18) with ε corresponding to a5% variation, computed for the parameters used in the numerical simulations and marked ascalibrated in Table 2. The average values have been taken on 100 independent runs, in thecase of oxygen at 21%. Observed variables N and the CSp diameter have been considered atthe final time of 72 h. In the first row for N and for the diameter the reference average valueand the standard deviation are shown.Parameters changed N SV N diameter SV diameter average=73.31 average=123.28 µ m(st. dev.=13.42 %) (st. dev.=5.71%) α + ε -0.68 % 0.14 -0.21% 0.04 α − ε -1.28 % 0.26 -0.83% 0.17¯ c ( ϕ ) + ε -2.08 % 0.41 -1.38% 0.27¯ c ( ϕ ) − ε -1.28 % 0.26 -0.83% 0.16¯ c ( ϕ ) + ε +0.31 % 0.06 +0.36% 0.07¯ c ( ϕ ) − ε +0.29 % 0.06 -0.51% 0.10 H + ε -0.29 % 0.06 -0.30% 0.06 H − ε +2.74 % 0.55 -0.73% 0.14 k + ε +0.79 % 0.16 +0.40% 0.08 k − ε -0.28 % 0.05 -2.57% 0.51 k + ε +0.87 % 0.17 -0.05% 0.01 k − ε -0.25 % 0.05 -1.19% 0.24 σ ( ϕ ) + ε -1.36 % 0.27 -1.15% 0.23 σ ( ϕ ) − ε -1.87 % 0.37 -1.43% 0.29 σ ( ϕ ) + ε -2.63 % 0.53 -0.83% 0.17 σ ( ϕ ) − ε +0.87 % 0.17 -0.83% 0.17 ρ + ε +0.80 % 0.16 -0.70% 0.13 ρ − ε -0.53 % 0.10 -0.76% 0.15 λ ( ϕ ) + ε +1.32 % 0.26 -1.55% 0.31 λ ( ϕ ) − ε -2.55 % 0.51 -0.84 % 0.16 λ ( ϕ ) + ε +1.80 % 0.36 -0.15% 0.03 λ ( ϕ ) − ε -2.01 % 0.40 -1.01 % 0.20 We performed numerical simulations of growth and differentiation of CSp considering the O and TGF- β as the only key-regulatory biological mechanism. We compared the compositionand the structure of the CSp at two typical experimental oxygen conditions.First, we analyzed the normal culture condition which corresponds to the 21%concentration of the oxygen O . The state of the CSp at the initial and three followingtimes is presented at Figure 5. Differentiation levels are marked by different colors: green forthe least differentiated cells ( ϕ = 1), blue color is for the intermediate level of differentiation( ϕ = 2), while the red color labels cells with the highest degree of maturation ( ϕ = 3). The16eported results are in good agreement with the biological observations. We reproduced theCSp biological system consisting of a central core of less differentiated but faster proliferatingcells surrounded by more specialized ones. The oxygen and TGF- β concentrations at differenttimes are shown at Figure 6.Then we considered the hypoxic culture conditions in the CSp environment, that is the5% oxygen concentration. The corresponding structure and composition of CSp, shown atFigure 7, match the real situation. The size of the sphere is larger than that experimentallyobserved, because of the equally circular shape of cells and their constant radius in our model.To compare and summarize the above results we show at Figure 8(a) in the form of piecharts the cell composition after 72 h in the two above cases. Similarly, Figure 8(b) presentsthe growth of CSp diameter, and Figure 8(c)–8(d) the time evolution of the ratio betweenthe total mass at time t and the initial mass, m rel ( t ) := m ( t ) /m (0), respectively for oxygenand TGF- β . We observe that the oxygen levels decrease more at the hypoxic conditions (5%)due to the increased proliferation and, as a result, initially the diameter of the sphere at 5%is larger than at 21%. However, the increasing deficiency of the oxygen lead to the formationof the necrotic core and the slower growth. Because of the initial condition equal for eachcell and equal maturation time, cells undergo a higher degree of synchronisation. As a resultwe observe the characteristic steps in the curves describing the growth of the CSp diameter,in correspondence of the beginning of the proliferation cycle.The proposed model appears to be in a good agreement with the observed biologicalexperiments at hand (Messina et al, 2004; Chimenti et al, 2012). We were able to reproducethe basic structures of CSps at two different oxygen concentration levels. In case of21% of the environmental concentration we observed the central region with proliferating,undifferentiated cells and the external ring of less multiplicating but differentiating cells. Onthe other hand, at 5% an evident necrotic core is observed, without many cells at higherdifferentiation level in the surrounding layers.17
50 100 150 200050100150200 (a) (b) t = 0 ,N = 15 ,N = 0 ,N = 0 ,N d = 0 t = 24 h,N = 5 ,N = 10 ,N = 0 ,N d = 0 (c) (d) t = 48 h,N = 9 ,N = 31 ,N = 4 ,N d = 0 t = 72 h,N = 19 ,N = 49 ,N = 15 ,N d = 0 Figure 5: (a)–(d) Numerical simulation of the CSp growth for a 21% oxygen environmentalconcentration, at times t = 0 , , ,
72 h. System (17) is solved in a domain Ω = [0 , × [0 , µ m ), with the model parameters given by Table 2. Differentiation levels are markedby different colors: green for the least differentiated cells ( ϕ = 1), blue color is for theintermediate level of differentiation ( ϕ = 2), while the red color labels cells with the highestdegree of maturation ( ϕ = 3). Values N , N , N , N d , in the subplot captions indicate thenumbers of cells in the CSp, respectively for each state of maturation and for dead cells, ϕ = 1 , , ,d . 18 a) (b) (c)(d) (e) (f) Figure 6: Spatial distribution of the oxygen (a)–(c) and the TGF- β (d)–(f) in the case of21% oxygen environmental concentration, at t = 24 , ,
72 h. Equations (17) , are solved ina domain Ω = [0 , × [0 , µ m ), with the model parameters given by Table 2.19
50 100 150 200050100150200 (a) (b) t = 0 ,N = 15 ,N = 0 ,N = 0 ,N d = 0 t = 24 h,N = 16 ,N = 7 ,N = 0 ,N d = 0 (c) (d) t = 48 ,N = 41 ,N = 16 ,N = 0 ,N d = 7 t = 72 h,N = 114 ,N = 2 ,N = 0 ,N d = 36 Figure 7: (a)–(d) Numerical simulation of the CSp growth for a 5% oxygen environmentalconcentration, at times t = 0 , , ,
72 h. System (17) is solved in a domain Ω = [0 , × [0 , µ m ), with the model parameters given by Table 2. Differentiation levels are markedby different colors: green for the least differentiated cells ( ϕ = 1), blue color is for theintermediate level of differentiation ( ϕ = 2), while the red color labels cells with the highestdegree of maturation ( ϕ = 3). Black region indicates a necrotic/apoptotic core composed bydead cells. As in Figure 5, N , N , N , N d , in the subplot captions show the composition ofthe CSp at the displayed times. 20 a) Time (h) C Spd i a m e t e r ( µ m ) O at 21% O at 5% (b) Time (h) R e l a t i v e m a ss o f O O at 21% O at 5% (c) Time (h) R e l a t i v e m a ss o f T G F - β O at 21% O at 5% (d) Figure 8: (a) Pie charts describing the relative percentage distribution of the three typesof cells (e.g. ( N /N tot ) ∗ β relativemasses, m rel ( t ) := m ( t ) /m (0) versus time. For all the figures blue and red curves correspondsrespectively to 21% and 5% oxygen concentration. In figures (b)–(d) the solid line representsthe average value computed over 100 simulations, while the dotted curves show the intervalcorresponding to one standard deviation. 21 Conclusions and perspectives
In this paper we developed a hybrid mathematical model describing the dynamics of cardiacbiopsy-derived stem cells leading to the formation of the so-called Cardiosphere from a clusterof a few cells.The complexity of the problem is extremely high and the biological phenomenaleading to the formation of a CSp are still not completely understood. Moreover thequantitative characterization of the processes misses the estimate of some parameters essentialfor the accuracy of a mathematical model. Nevertheless, the proposed mathematical modelappears in perfect agreement with the expected biological and pharmacological applicationshown in Messina et al (2004); Chimenti et al (2012, 2010); Forte et al (2012). In particular,we were able to reproduce the basic structures of CSps assuming that only two elementsare sufficient to describe the overall growth and maturation of the CSps: the oxygen asnutrient’s emblem and the TGF- β as the chemical differentiation signal representative. Ourchoice to select these two key regulators of the CSps growth and differentiation gives us manymore chances to easily foresee the main possible cellular consequences following experimentalbiological or pharmacological modifications, with a strong impact in both research qualityand cost. The layered growth of the CSps, with a central proliferating region, surroundedby differentiated and less proliferating cells, has been observed in the numerical experiments.In this regard, the biological technology has a high potential for therapeutic purposes andpossible strategies for further differentiated cells are under investigation (Lee et al, 2013; Chenet al, 2014). Furthermore, in a possibly more general framework, our modeling tool couldbe applied to many other cellular spheroid culture systems, such as tumor cells spheroids,induced-progenitor cells, embryoid bodies-derived stem cells. References
Agathocleous M, Harris WA (2013) Metabolism in physiological cell proliferation anddifferentiation. Trends in Cell Biology 23(10)AIS Project (2001) Diffusion of molecules. Web page, Department of Mathematics, Universityof British Columbia, URL:
Anderson ARA, Chaplain MAJ, Rejniak KA (eds) (2007) Single-cell-based models in biologyand medicine. Mathematics and Biosciences in Interaction, Birkh¨auser, Basel, SwitzerlandAnderson DF, Ermentrout B, Thomas PJ (2015) Stochastic representations of ion channelkinetics and exact stochastic simulation of neuronal dynamics. Journal of ComputationalNeuroscience 38(1):67–82Bayly PV, Taber LA, Carlsson AE (2012) Damped and persistent oscillations in a simplemodel of cell crawling. J R Soc Interface 9(71):1241–1253Bell GI, Dembo M, Bongrand P (1984) Cell adhesion. Competition between nonspecificrepulsion and specific bonding. Biophys J 45(6):1051–1064Brown G, Hughes PJ, Michell RH (2003) Cell differentiation and proliferation: simultaneousbut independent? Exp Cell Res 291:282–28822haplain MAJ (2011) Multiscale mathematical modelling in biology and medicine. IMAJournal of Applied Mathematics 76(3):371–388Chen WP, Liu YH, Ho YJ, Wu SM (2014) Pharmacological inhibition of tgf receptor improvesnkx2.5 cardiomyoblast-mediated regeneration. Cardiovascular Research 105(1):44–54, DOI10.1093/cvr/cvu229Chimenti I, Ruckdeschel Smith R, et al (2010) Relative Roles of Direct Regeneration VersusParacrine Effects of Human Cardiosphere-Derived Cells Transplanted Into Infarcted Mice.Circul Res 106:971–980Chimenti I, Gaetani R, Barile L, Forte E, Ionta V, Angelini F, Frati G, Messina E, GiacomelloA (2012) Isolation and Expansion of Adult Cardiac Stem/Progenitor Cells in the Form ofCardiospheres from Human Cardiac Biopsies and Murine Hearts. Met Mol Biol 879:327–338, DOI 10.1007/978-1-61779-815-3 19Clarelli F, Di Russo C, Natalini R, Ribot M (2015) A fluid dynamics multidimensional modelof biofilm growth: stability, influence of environment and sensitivity. Math Med Biol DOIdqv024v1-dqv024Di Costanzo E, Natalini R (2015) A hybrid mathematical model of collective motion underalignment and chemotaxis. Preprint, ArXive Repository, URL: http://arxiv.org/abs/1507.02980
Di Costanzo E, Natalini R, Preziosi L (2015) A hybrid mathematical model for self-organizingcell migration in the zebrafish lateral line. J of Math Biol 71:171–214D’Orsogna MR, Chuang YL, Bertozzi AL, Chayes LS (2006) Self-Propelled Particles withSoft-Core Interactions: Patterns, Stability, and Collapse. Phys Rev Lett 96(10):104302Enderling H, Hlatky L, Hahnfeldt P (2012) The promoting role of a tumour–secretedchemorepellent in self-metastatic tumour progression. Mathematical Medicine and Biology29:21–29Fletcher AG, Murray PJ, Maini PK (2015) Multiscale modelling of intestinal cryptorganization and carcinogenesis. Math Mod Meth Appl S 25:2563–2585Forte E, Miraldi F, Chimenti I, et al (2012) TGFb-Dependent Epithelial-to-MesenchymalTransition Is Required to Generate Cardiospheres from Human Adult Heart Biopsies.Stem cells and Development 21(17)Fournier MF, Sauser R, Ambrosi D, Meister JJ, Verkhovsky AB (2010) Force transmissionin migrating cells. J Cell Biol 188(2):287–297Itoh M, Chitnis AB (2001) Expression of proneural and neurogenic genes in the zebrafishlateral line primordium correlates with selection of hair cell fate in neuromasts. MechDevelop 102:263–266Joie J, Lei Y, Colin T, Durrieu MC, Poignard C, Saut O (2013) Modelling of migration andorientation of endothelial cells on micropatterned polymers. Research Report RR-8252,Inria Institute, URL: https://hal.inria.fr/hal-00795238