CC BY
14Chemo-Mechanical Elastic Modeling of Carcinoma Development
D.A. Bratsun1*, A.P. Zakharov2, and L.M. Pismen2
1 Perm National Research Polytechnical University, Department of Applied Physics, Perm, 614990 Russia 2 Faculty of Chemical Engineering, Technion—Israel Institute of Technology, Haifa 32000, Israel
We propose a multiscale chemo-mechanical model of the cancer tumor development in the epithelial tissue. The epithelium is represented by an elastic 2D array of polygonal cells with its own gene regulation dynamics. The cell morphology was calibrated in the model based on the experimental data. The model allows for the simulation of the evolution of multiple cells interacting via the chemical signaling or mechanically induced strain. The used algorithm takes into account the division and intercalation of cells, as well as the transition of normal cells into a cancerous state triggered by a local failure of spatial synchronization of cellular rhythms driven by transcription/translation processes. Both deterministic and stochastic descriptions of the system are given for chemical signaling. The simulations reproduce a distinct behavior of invasive and localized carcinoma. Generally, the model is designed in such a way that it can be readily adjusted to take into account any newly understood gene regulation processes or feedback mechanisms affecting chemo-mechanical properties of cells.
Keywords: cancer modeling, signaling, circadian rhythms, time-delay, systems biology
1. Introduction
Mathematical modeling of cancer has been growing immensely as one of the challenges of mathematics and physics applied to biology and biochemistry. The principal difficulty in modeling of cancer (as, in fact, of any biological system) is the multiscale nature of the phenomenon [1]. One can identify at least three natural scales, with different stages of the disease development. Processes on the cellular scale are triggered by signals stemming from the subcellular level and have an impact on the macroscopic scale, i.e. on the organism as a whole, when tumor grows and spreads [2].
At the cellular scale, a system of coupled ODEs (Ordinary Differential Equations) can be used to model large cell populations, where each variable corresponds to a well-defined biological property characteristic of all cells of the same population. This approach has been further developed to account for more fine effects (see, for example, review paper [3]). The advantage of this approach is that the models are easily tractable, and enable a relatively fast identification of the parameters. On the downside, however, the
* Corresponding author
Prof. D.A. Bratsun, e-mail: [email protected]
models do not account for potentially important phenomena, such as spatial aspects and heterogeneity among cells. Another type of modelling, supplementing the population dynamics with additional variables determining the average structure of a population of cancer cells (for example, the age of cells), can be assihned to the class of semi-phe-nomenological models with internal structure [4].
A large body of literature has been devoted to models linking the cellular scale to the macroscopic tissue scale [2]. Most commonly, the tumor is considered as a continuous medium.
It is described by a system of partial differential equations including the mass balance equation for the cellular medium and reaction-diffusion equations describing the field of chemical signals exchanged between cells. For example, the tumor may be viewed as a porous matrix [5] which interacts with a filling liquid (healthy cells). In a very recent work of this kind [6], the tumor was represented as a multiphase medium simulated in a 3D geometry by the finite elements method. The disadvantage of such models is their inherently phenomenological character.
An alternative approach is cell-based discrete modeling. Unlike continuum models, discrete models have the ability to track the behaviour of individual cells. Due to advances in biotechnologies, there is an increasing amount
Original text © D.A. Bratsun, A.P. Zakharov, L.M. Pismen, 2016
© Institute of Strength Physics and Materials Science, Siberian Branch, Russian Academy of Sciences, 2016. All rights reserved.
of experimental data available at a single cell level that can be used for improving the mathematical models. Early models including single cell data were based on cellular automata or random walk of cells. A large number of models combining interactions at sub-cellular and cellular level have been reported [2]. The small number of reactant molecules involved in gene regulation can lead to significant fluctuations in intracellular mRNA and protein concentrations and the consequences of such noise at the regulatory level have been analyzed in numerous recent studies [7]. Considerable experimental evidence exists that stochasticity plays major role in gene regulation, both due to intrinsic and extrinsic factors. There are several different approaches to modeling of stochastic chemical reactions: direct Gillespie simulations [8], exact master equation analysis and simplified descriptions based on the Langevin equation [7]. In addition, the transcriptional and translational processes are known to be compound multistage reactions involving the sequential assembly of long molecules. This can provoke a time lag in gene regulation processes. The combined effect of time delay and intrinsic noise on the temporal dynamics was explored in [9, 10].
The multiscale nature of cancer requires mathematical modeling approaches that encompass different biological scales. Due to the complexity of such approaches the number of literature reports is limited. One of the first attempts to implement a hybrid approach included a model of cellular automata whose state is determined by a continuous distribution of oxygen around a blood vessel near the origin of the tumor [11]. In another remarkable work [12], a spherically growing tumor within the framework of a lattice model of mechanically interacting discrete cells was examined. A dynamic system describing the processes of gene regulation has been applied separately to each cell in the population. Another multiscale chemo-mechanical model of cancer tumor development in an epithelial tissue has been re-
cently proposed in [13, 14]. The model is based on transition of normal cells into the cancerous state triggered by a local failure of spatial synchronization of the circadian rhythm. A good review of recent developments in multiscale cancer modeling can be found in [1].
In the present paper we discuss the main principles of the tumor growth modeling within a cell-based mechanical model of deformations and rearrangement of epithelial tissue coupled with the transcriptional regulations and intercellular signaling earlier developed in [13-9]. The main hypothesis is that tumor formation is driven by disruption of the circadian rhythm in the epithelial tissue. Circadian rhythms are common to almost all living organisms. The timing of circadian clocks is established in a cell-autonomous manner by a self-sustaining molecular oscillator that includes intertwined negative and positive transcription/ translation-based feedback loops. It has been recognized in recent years that core circadian genes are important in tissue homeostasis and tumorigenesis. Many studies have shown (see, for example, [20, 21]) that disruption of the circadian clock is implicated in gene deregulation leading to the development of cancer and other diseases.
2. Model description
2.1. General Principles
Epithelium can be defined as a relatively avascular aggregation of cells that are in apposition over a large part of their surfaces, and are specialized for absorptive, secretory, protective, or sensory activities. The cells of lining and covering epithelium (as in the lining of intestine) are arranged in sheets whereas glandular epithelium consists of complex aggregates of epithelium cells. We focus on the first type, and construct a 2D model of a single layer of epithelial cells. Different shapes can be distinguished in microscopic sections, but the distinctions are often blurred.
Fig. 1. (a) - Photomicrograph of the papillary thyroid carcinoma demonstrating nuclear overlapping, (b) - Photomicrograph of the Pap test smear showing the presence of cervical cancer cell (a larger nucleus surrounded by normal cells). Both micrographs were taken in the Perm State Clinical Hospital.
The cells adhere to each other forming specialized attachment structures (desmosomes) that ensure coherence and strength of tissue. We will further concentrate on the modelling of carcinomas. Carcinoma is a type of cancer that develops from epithelial cells when DNA is altered or damaged to such an extent that the cells start to exhibit abnormal malignant properties. It may affect skin, breast, lung, prostate, and colon, and is among the most common types of cancer in adults.
Figure 1 shows two different types of carcinoma cells found in patients of Perm State Clinical Hospital. The first is a papillary thyroid carcinoma which develops from epithelial cells in the thyroid gland (also known as follicular cells), responsible for the production and secretion of thyroid hormones. Follicular cells are of the cuboidal type, and are arranged in spherical follicles (Fig.1a). It can be seen that large, crowded, overlapping and sometimes empty-looking nuclei cells are the characteristic features of this carcinoma. Another type of carcinoma cells is shown in Fig. 1 (b) presenting the Pap smear of a cervical cancer patient. The invasive cervical carcinoma cells are mostly squamous; they are relatively small in size and have large nuclei. This is in contrast to normal epithelial cells (also shown in Fig. 1 (b)) that are larger, almost rectangular in shape and have smaller nuclei.
The following key features make our model suitable for the realistic simulations of the epithelium:
- Cells change size and shape in the process of tissue evolution.
- Tissue spreads by the mechanism of cell division.
- Individual cells move within the tissue by the mechanism of intercalation.
- Neighboring epithelial cells exchange chemical signals through their common borders.
- Dynamics of signaling species take part in the regulation of intracellular processes.
- Normal cells are able to transform into cancer cells.
- A new species of cells with its own set of physical and mechanical properties.
The following properties of cancer cells are taken into account:
- Cancer cell never undergo the reverse transition to normal cells.
- Cancer cell exhibits a number of alterations on cell surface in the cytoplasm and in their genes.
- Cancer cell show uncontrolled mitotic divisions causing disorganized growth.
- A tumor can be formed due to uncontrolled growth and division of cancer cells.
- Cancer cells are far less adhesive than the normal cells, and therefore tend to wander through tissues causing cancerous growth in different parts of the body.
In order to obtain a realistic description of the epithelium that reproduces the irregular stressand velocity distribution, we developed a cell-based model combining ba-
Fig. 2. Principal components of the model and their interrelation.
sic mechanical and reaction-diffusion processes. The main components of our model are mechanical, involving elastic interactions between cells and their spreading, and genetic, involving cell transformation and signalling influenced by circadian rhythms. Their relationships are schematically presented in Fig. 2. Below we discuss each element of the model in more detail.
2.2. Cell morphology
In order to calibrate the mathematical model, cytology materials obtained from patients of Perm State Clinical Hospital were examined. Distribution histograms of physical parameters reflecting cancer and normal cell morphology— characteristic size, perimeter and area were obtained from a large number of cell micrographs (see a representative example in Fig. 1b). Figure 3 shows the distribution histograms of normal (green) and cancerous (red) cells on the perimeter. It can be seen that the cancer cells have a smaller average size. This can be explained by an increased internal pressure in a tumor that arises due to rapid cell division.
It is worth noting that our results obtained by direct measurements of morphological parameters from photomicrographs are in good qualitative agreement with cell morphology measurements performed using new sophisticated optical techniques [22].
2.3. Chemo-mechanical cell-based model of epithelial tissue
The mechanical model presents the epithelium as an elastic two-dimensional array of cells, approximated by polygons. We take as the initial configuration a regular hexago-
Fig. 3. Distribution histograms of the normal and cancer cells on the perimeter.
nal lattice. In the course of spreading, division and transformation, it becomes distorted and incorporates also polygons with a different number of vertices. Controlled spreading and deformation of epithelial layers, two-dimensional sheets of cells that tightly adhere to one another, is a key process in both adult and developing tissues. During embryogenesis, spreading and folding of epithelia plays a central role in gastrulation, an event that initiates the formation of three-dimensional structures of tissues and organs, as well as in the subsequent dorsal closure. Mechanisms of collective cell migration are a subject of intensive research [23]. The process is mediated by the intracellular actin machinery and proteins of adherent junctions, and is regulated by complicated chemical pathways, that still remain largely unexplored, and are influenced both by genetic and extrinsic factors.
Epithelial tissue is a layer of cells covering the surface of an organ or body. Thus, we use only quasi-two-dimensional system in modeling the epithelium behavior, that makes the calculation easier [13-16]. The cells always remain attached to each other forming a continuous two- dimensional epithelial surface. The curvature of the layer is presumed here to be small compared to the cell size and can be neglected (a more complicated case of the deformed layers of cells has been considered recently in [18]).
The mechanical model is based on the elastic potential energy U of the tissue:
u=2 E (¡P2 +n(A - A0)2), (1)
2 cells
where P and A stand for the perimeter and area of cell respectively. Here the coefficient ¡x characterises the effect of contractile forces within the perimeter of the cell, n characterizes the elastic resistance to stretching or compressing the cell with respect to the reference cell area A. The vertices of the polygons representing cells form a lattice. Evolution of the tissue occurs by moving the lattice nodes. We
define the mechanical force acting on any jh node as
=-f ■ ®
where Rj denote the position of the node.
We consider the internal movement of cells in the tissue as a strongly overdamped process, so the equation of motion for i-th node can be written as
V = KFtH (IFil-Fo),
at
(3)
where Vi is the velocity, K is the mobility coefficient, His the Heaviside function. The threshold force F0 has been introduced in (3) to take into account the situation when the node remains immobile even if force is not zero.
An important feature of the model is the ability of cells to divide. This allows the tissue to carry out internal movement through the redistribution of internal stresses in the environment. It is assumed that the division occurs when the longest edge of the polygon and its corresponding opposite edge are divided in half (Fig. 4a). In order to minimize the growing disorder in the distribution of nodes in the process of cell division, the probability to divide was connected to the number of vertices of the polygon n according to the following formula [16]:
Pdiv = P0<T6, (4)
where q is a distribution constant and P, is a scaling factor. One can see that for q > 1, the polygons with a number of sides over 6 experience the division more frequently. Thus, the most likely shape of the cell according to (4) is a hexagon. For cancer cells, the dependence of the number of nodes is suspended and the division rate is set at considerably higher level.
Another mechanism which increases the liquidity of the tissue and excludes the severe deformation of some cells is intercalation [13, 16]. We introduce here a special parameter l0 which determines the moment when the intercalation can occur. Then the probability of the event can be written in the simple form:
J1, li < ¿0>
[0, li > l0.
One can see from (5) that if the length of the border separating two cells becomes less than l0, it is substituted by a link of a slightly larger length in the normal direction (Fig. 4b). The intercalation is known to be important in many
P =
1 int
(5)
Fig. 4. Elements of chemo-mechanical model of epithelium: cell division (a); cell intercalation (b).
tissue reshaping processes. Since cancer cells are far less adhesive than normal cells, the respective l0 is set at a higher level, thereby making the intercalation of cancer cells more probable. Altogether, Eqs. (1)-(5) define the mechanical dynamics of tissue on the cellular level.
2.4. Genetic model of circadian rhythms: single-gene auto-repressor model with dimerization
The key feature of our genetic model is the influence of periodic inputs. Biological rhythms are periodically repeated changes, which are quite characteristic for living matter on every level of its hierarchy, starting from the molecular and subcellular levels and up to the biosphere as the whole. In a living organism, rhythms are closely connected with their adaptation to the environment during the evolution. We concentrate here on circadian rhythms that synchronize with daily changes of the environment. A remarkable feature of these rhythms is that they are not simply a response to the 24-hour environmental cycle imposed by the Earth's rotation, but are generated internally by cell-autonomous biological clocks [24]. After decades of research, the genetic mechanism of circadian oscillations has been widely recognised as a core of this phenomenon. It was realized that their mechanism works even on the scale of one or several genes by revealing itself in RNA and protein fluctuations in transcription and translation processes. As soon as transport proteins pass through the cell membranes and trigger intracellular interactions, circadian oscillations inevitably develop at the intercellular scale. At the organism scale, the signals from separate cells should be synchronised, thus developing unified rhythms for the whole organism. For example, in mammals, it is believed that a master pacemaker in the hypothalamus orchestrates temporal alignment of behaviour and physiology by transmitting daily signals to multiple clocks in peripheral tissues. The molecular mechanisms of circadian rhythms for some organisms are already well understood [25], but the studies
largely concentrate on the temporal rather than spatial organization of rhythms.
Since there is extensive experimental evidence that gene deregulation influenced by circadian clock is implicated in the development of cancer [20, 21], the molecular mechanism of the circadian rhythm disruption and its propagation due to cell-to-cell signaling should become an essential factor in transformation of cells (Fig. 2). We assume that a normal cell can turn into a cancer cell due to a local disruption of the circadian rhythm in the peripheral epithelial tissue. For this reason, the mathematical models of rapid synchronization of the entire community of oscillators developed for SCN cells (see, for example, [26]), are not very suitable for our purposes. Instead we apply a single-gene auto-repressor model with dimerization where the negative feedback loop is delayed in time as suggested in [9]. Time-delay seems to be the most common cause of oscillations in genetic systems, since gene regulation processes are typically very slow and comprise multistage biochemical reactions engaging the sequential assembly of long molecules, and therefore are likely to generate time delays.
Let us consider a single-gene protein synthesis with negative auto-regulation (Fig. 5a). This is a popular motif in genetic regulatory circuits, and its temporal dynamics has been analyzed within both the deterministic and stochastic framework [7]. The generalized version of this system taking into account that the production of the auto-repressor protein takes a finite amount of the delay time has been studied in [9, 10].
Suppose that protein can exist both in the form of monomers Xand dimers Y. The transitions between them with the
rates k+ d and k-d are X + X ——— Y, Y-
->X + X.
(6)
We assume also that the protein can be degraded with the rate B and produced with the rate A respectively:
X —^0, DO
Di+T+ Xt
(7)
200
240 260 280 Time, h
320
Fig. 5. Network architecture of the molecular components of the circadian rhythm: a single-gene auto-repressor model with time-delay where only x gene defines the oscillatory activity (a). X stands for the concentration of protein produced by x gene; a typical time series obtained during single-cell simulation within deterministic description given by (9): the dynamics of X protein with parameters t = 8 h, A = 5000 nM/h, B = 5.0 1/h, 8 = 0.1 1/nM, Ô = 0.2 1/nM (b)
The synthesis occurs at time t + t if the chemical state of the promoter site of the x gene at time t is unoccupied (Do). Otherwise (Dl), the production is blocked. The transitions between the states occur with rates k+1 and k_1 during binding and unbinding of some dimer are respectively D0 + Y -1 > D1, D1 - > D0 + Y. (8)
In order to describe the intercellular signaling, we assume that the X monomers can be transported via cells membranes. Thus, Eqs. (6)-(8) define the kinetics of gene regulation both the single cell level and whole tissue level.
2.5. Genetic model of circadian rhythms: deterministic description
The main approximation is that the reactions of dimer-ization (6) and binding/ unbinding (8) are fast in comparison with production/degradation of protein (7). Thus, its dynamics has to enter quickly into a local equilibrium and we obtain
5 a 6
dxi (t ) dt
aN
1
A
1 + esx2(t _ t)
1 + 4ex (t )
E Pj(Xj (t) _ X (t)).
jGadj(i)
_ BXi (t )
(9)
where the subscripts refer to cells, x stands for the concentration of the free monomers of the X protein, a is the transfer coefficient, N is the copy number, 8 = k+ d/k-d , £ = k+1/k-1 and adj(i) stands for "adjacent to ith-cell". We assume for the simplicity that the same X protein transports the circadian signal outside the cell. It is assumed that the monomers of the Xprotein are transported diffusively from one cell to the other, whereas its flux does not depend on the distance between the two cells i and j but is proportional to the boundary length Pj. This implies that the transport is limited by the transfer though cell membranes. The link between sub-cellular and macroscopic scales is established through the Eq. (9), since the x field is global for the whole tissue. After a cell divides, the daughter cells inherit the phase of the circadian rhythm of the parent cell.
The neutral curve for the Hopf bifurcation of (9) within the deterministic approach at a single- cell level was derived in [9]. The numerical study reveals the oscillatory dynamics above the bifurcation as it is shown in Fig. 6b.
In order to study the spatial effects within the deterministic description, the set of delay differential equations (9) has been solved using the explicit Euler method, whose stability was warranted by a sufficiently small time step. This procedure was synchronized with the simulation of the mechanical evolution governed by Eqs. (1)-(5). The initial configuration of the system is a hexagonal lattice comprising 1560 cells with random phase distribution. The tissue as a whole has the form of a stripe. Figure 6a presents the results of numerical simulation of the X protein pattern with parameters taken above the Hopf bifurcation. The nonlinear dynamics includes the slow development of spiral
traveling wave pattern which arises against the synchronized field oscillating in the background. The oscillation period is approximately equal to the triple delay time.
2.6. Genetic model of circadian rhythms: stochastic description
The small number of reactant molecules involved in the gene regulation can lead to significant fluctuations in protein concentrations. Since the pioneering works published in the early 2000s, there have been numerous studies dealing with the influence of such noise at the regulatory level [7]. A review of recent developments in the field can be found in [27]. Additional point to emphasize is that the transcription-translation processes are compound multistage reactions involving the sequential assembly of long molecules. This may result in a time lag in gene regulation processes. When the delays are short compared with other characteristic time scales of the system, one can safely ignore them in simulations. However, if the lags become longer than other processes, the system has to be considered as non-Markovian, and one should account for this in both deterministic and stochastic descriptions. The joint effect of the intrinsic noise and time delay on the temporal behaviour during gene regulation have been studied first in [9, 10]. In [9] we have suggested a generalization of the Gillespie algorithm [8] widely used to simulate statistically correct trajectories of the state of a chemical reaction network that accounts for delay. Based on this technique, we have shown that quasi-regular fluctuations can arise in the stochastic system with delay even when its deterministic counterpart exhibits no oscillations [10]. Lately a large number of works have been published developing this research line (see a recent review [28]), mainly focusing on further improving the algorithm and on studying the temporal dynamics of gene systems with delays.
In [29] it was stated that "space is the final frontier in stochastic simulations of biological systems". The problem is that despite the considerable body of spatio-temporal experimental data has been accumulated, the stochastic models of biochemical processes focus mostly on temporal dynamics. If in the past years considerable progress has been made in spatial stochastic simulations of Markovian processes [27, 29], however the studies of non-Markovian stochastic systems are still very rare. The theoretical difficulties seem to be clear: the generalization of the Gillespie algorithm to the case of the spatial dynamics of time-delayed processes is still waiting for solution. In order to describe the spatial stochastic effects, here we use a hybrid model, which is constructed as follows. The dynamics of the proteins in each cell has been obtained by performing direct stochastic simulations of the reactions (6-8) using the modified version of the Gillespie algorithm [9]. The signaling between cells is still organized as diffusive transport of the X monomers from one cell to the other according to finite-difference formula:
Xt+At = Xt + G (Xt, Xt) +
At £ aLiJ( Xj, Xt) , (10)
jGadj (i)
where [...] stands for the integer part of the expression and G denotes the generalized Gillespie algorithm. The time step in (10) is equal to the time step for the integration of the mechanics ofthe system (1)-(5): At = 0.05. Since the typical time step of the stochastic system is much less, one needs to dock the numerical schemes for the tissue mechanics and stochastic dynamics of the protein.
Figure 6b presents the stochastic pattern formed by the Xmonomers in the tissue with the same parameter values as in Fig. 6a. We found that nonlinear dynamics of spatially extended system consists of two distinct oscillatory modes, just like it was in the deterministic case. One is a quasistanding wave pattern oscillating with the period equal to the triple delay time. The second oscillatory mode consists of traveling waves which arise from selected initial disturbances. In fact, the stochastic pattern looks very similar to its deterministic counterpart obtained for the same parameters (compare the frames in Fig. 6). The wavelength of the structure is found to depend on the copy number N of signaling whose growth enhances fluctuations and diffusive fluxes between cells.
2. 7. Transformation of cells
As mentioned above, the main circadian genes appear to strongly influence tumorogenesis. The rhythmicity demonstrated in the expression of clock-controlled genes regulates various functions of cells, including their division and proliferation. Resynchronization of this rhythmicity can be involved in some pathologies, including the development of tumours. There is increasing evidence linking malfunction of the bioclock work with pathogenesis of cancer [20,
21]. In a review paper [30], a large number of examples ofthe connection between circadian genes, circadian periodicity, aging-related phenotypes, and cancer are given. An important result obtained experimentally was reported in [31]. The authors have checked whether the circadian rhythm was connected to the cell-cycle oscillation in immortalized rat-1 fibroblasts by observing cell-cycle gene promoter-driven luciferase activity. It was revealed that there was no direct phase relationship between the circadian and cell rhythms. These data imply that the circadian system does not govern the cell-mitosis rhythm in rat-1 fibroblasts. Thus, it was suggested that there is no direct coupling between the circadian rhythm and cell cycle but the timing of cell mitosis is synchronized with the rhythmic host environment. Based on these and other studies, the main idea of the alteration mechanism is a synchronization failure of a local oscillation phase in the common field of spatially synchronised circadian rhythms in the epithelial tissue.
Preliminary numerical simulations of the spatially extended problem (1 )-(5) with the circadian model (6)-(8) have shown that when the number of cells is large enough, a complete synchronisation, i.e. the total alignment of the oscillation phases in all cells, cannot be achieved when the coupling is small. Instead, the cells are organised in collective spatio-temporal patterns including clusters of cells oscillating with almost the same phase. The thin layer of cells between clusters exhibit oscillations with intermediate phase. Thus, we introduce a new variable which characterises the local dephasing defined as the phase difference of the cir-cadian rhythm of an ith cell and the average phase of the adjacent cells:
(t ) = ^
X (t ) Xk (t )
max( Xt ) max(Xk )
(11)
kGadj (i )
where phase values are normalized by its maximal value. If
Fig. 6. X protein pattern in the epithelium formed of 1560 cells within (a) deterministic and (b) stochastic description, with parameters A = 5000 nM/h, B = 5.0 1/h, ô = 8 h, N= 1, < = 0.05 L/h, k+d = 200 1/nMh, k-d = 1000 1/h, k+i = 100 1/nMh, k-i=1000 1/h. Both frames correspond to time 340 h
the ith cell is in a fully synchronised field, the value of Oi, obviously, will be zero. On the other hand, the dephasing (11) increases near the boundary of the clusters. The maximum value of the dephasing would occur in case of an isolated cell within a phase cluster. According to our hypothesis, this cell is most at risk to become a cancer cell. Apparently, in a dynamically changing pattern where cluster boundaries are varying in time, the cell transformation should occur much less frequently.
Then we can introduce the phenomenological state equation for the ith cell [13, 16]: dZi (t)
dt
■ = —Z (t )(1 - Z (t ))(Z - Z (t )) + aOi- V (t ), (12)
where Z is the state function, X is the damping parameter, £(t) is an uncorrelated zero mean noise input in the range [-1; 1], a is the amplitude of the noise, assumed to be multiplicative. Without noise, Eq. (12) has two stable stationary solutions: Z = 0, standing for the normal state of the cell, and Z = 1 corresponding to the cancer state. These solutions are separated by an unstable steady state which separates their domains of attraction.
Thus, Eq. (12) constitutes the simple bistability model allowing the cells to transform into cancer state with some probability depending on the circadian pattern formation expressing the collective behavior of cells.
3. Numerical results
3.1. State function dynamics
The initial configuration of the system is set to be a regular hexagonal lattice comprising 1560 healthy cells and no cancer cells. The shape and location of each cell is defined by its nodes. The tissue as a whole has the form of a stripe with two free borders with periodic boundary conditions applied there. The typical values of the mechanical parameters for the normal cells areas follows: x = 1.0 F/L, n = 1.0F/l}, A = 3V3/2L2, K = 1.0L/(hF), F0 = 0.02F, P0 = 0.0002, q = 1.4, l0 = 0.15, X = 101/ h, Z = 0.15, a = 0.11/ (hL). Here spatial dimensions and force are measured in terms of arbitrary units L and F respectively.
The set of differential equations describing the dynamics of mechanical subsystem (1 )-(5), the circadian rhythms and signaling (6)-(8) (depending on the description (9) or (10)) and the state function dynamics (11)—(12) has been solved using the explicit Euler method, whose stability was warranted by a sufficiently small time step At = 0.005. The time step for the calculation of the molecular processes in cells was synchronised with the step of calculating the mechanical movement of the tissue cells. We take as the initial condition a random phase distribution. In the process of evolution, the rhythms in cells try to synchronise through a weak nonlinear interaction, generating the various spatio-temporal patterns. From the perspective of biology, an arbitrary distribution of phases in cells taken as the initial conditions looks artificial, since synchronisation of
Fig. 7. X protein concentration (a), dephasing O (b) and state function Z (c) in the epithelium composed of 1560 cells calculated within deterministic description at t = 350 h. System evolution starts from random phase distribution.
rhythms happens at the stage of embryonic development. Nevertheless, the numerical study with random initial conditions allows us to better understand the self- organisation properties of the system and to evaluate the mechanisms of pattern formation.
The typical values of the parameters governing the cir-cadian rhythms in each cell are as shown in Fig. 6. The next Figure 7 presents the results of numerical simulation of the temporal evolution of the system within the deterministic description. The patterns of the concentration x, the dephasing and the state function Z are shown for the time moment 350 h. We have found that the spatial dynamics includes the slow development of quasi-standing waves which arise against the mean field oscillating with the basic period 27 h (Fig. 7a). The dephasing calculated according to (11), (12) reaches a maximum value on the boundary of the standing waves domain (Fig. 7b). The cells that fall into this region, are at high risk to transform into cancerous state. Figure 7c confirms this conclusion: most of the transformed cells are grouped in a zone of uncertain oscillation phase characterized by the maximal dephasing. The important point is the small mobility of the pattern shown in the figure. Since the boundary of the standing waves change slowly, the same cells are always at risk of the transformation. Note that the parameter values for the transformation model (12) have been calibrated so that the number of transformed cells was sufficient to demonstrate the effect.
3.2. Tumor development
Figure 8 gives the example of the numerical simulation of a invasive tumour. The parametersdefining the properties of cancer cells are A = 3>/3/2 L2, ^ = 1.2 F/L, n = = 1.0 f/l3, l0 = 0.4. In order to focus upon the development of the tumour, the alteration process has been stopped during the simulation run just after the first cell turned into a cancer cell. This moment of time was fixed at t = 0 h. One can see that the size and shape of the cells becoming cancerous differ from those of the normal members of the community. Cancer cells are approximately twice as small and irregular in shape. The healthy cells that border on cancer cells experience a significant stress: they are squeezed and
I
shrunk under the onslaught of cancer cells. Since the period of the division of cancer cells is shorter, the tumour evolves rapidly, increasing the occupied area. The forced proliferation of cancer cells is the reason why in the bulk of tumour the cells are also strongly squeezed and may become irregularly shaped. One can notice that the threshold of intercalation l0 for cancer cells has been sharply increased in comparison with the normal cells. This condition allows to reduce a high level of potential energy in cells by restructuring their form due to the intercalation process. The characteristic rate of cell intercalation exceeds here the division frequency, and the front between the tumour and the healthy tissue becomes unstable. Because of the easy intercalation, active proliferation and relatively small size, the cancer cells actively change their location, move apart and migrate through the healthy tissue. Thus, this simulation models a key property of the invasive type of cancer being capable to penetrate into various tissues and organs. The instability of the front between normal and cancer cells, clearly visible in Fig. 8, looks similar to a fingering instability at the interface of two immiscible fluids when a fluid with a lower viscosity is pushed into a fluid of higher viscosity. Although tissue is neither a liquid nor a granular medium, and cannot be characterized as either miscible or immiscible, the analogy is natural, as the effective viscosity of the tumor is lowered by easy intercalation, while the driving pressure is generated by division of cancerous cells.
4. Conclusions
Cancer formation is a complex biophysical process, and its modelling requires a multiscale mathematical approach. In present review a minimal multiscale chemo-mechanical model which includes three natural scales of the tumor formation was proposed. Processes on the mesoscopic (cellular) scale are activated by circadian rhythm signals generated on the microscopic (subcellular) scale, and exert an
-0.2
-0.4
- -0.6
-0.8
Fig. 8. Numerical simulation of an invasive cancer tumor development in epithelial tissue: (a) - relative cell size with respect to the average size of a normal cell; (b) - state function Z where Z=1 and Z=0 stand for cancer and normal cells, respectively. The frames correspond to time 80 h.
influence over cooperative phenomena occurring at the macroscopic level characterised by a long-range coordination of cells. The medium where the tumour grows and spreads is modelled taking into account both transport of signaling species and mechanical interactions between the cells. In addition to detailed description of the model numerical results including simulation of tumor development are presented. It is worth noting that many mechanisms underlying tumor development currently remain unknown. Still, our model is designed in such a way that any newly discovered molecular mechanism can be easily integrated into it by using simple phenomenological relations. The model can be generalized to the case when epithelial cells cover a curved surface. Such situation is the most realistic one, because of the complex topology of epithelial tissue. Finally, the model can also be generalized to an arbitrary three-dimensional tissue. The number of epithelial cells that can simultaneously participate in tissue evolution is only limited only by the capacity of the computer.
Acknowledgments
The research has been supported by the Ministry of Education of Perm Region (grant C-26/004.4) and the Russian Fund for Basic Research (grant 14-01-96022r_ural_a).
References
1. Deisboeck TS, Wang Zh, Macklin P, Cristini V. Multiscale Cancer Modeling. Annu Rev Biomed Eng. 2011; 13: 127-155.
2. Bellomo N, Li NK, Maini PK. On the foundations of cancer modelling: selected topics, speculations, and perspectives. Math Mod Meth Appl S. 2008; 18: 593-646.
3. Kimmel M, Lachowicz M, Swierniak A. Cancer growth and progression, mathematical problems and computer simulations. Int J Appl Math Comput Sci. 2003; 13: 279-429.
4. Dyson J, Villella-Bressan R, Webb GM. The steady state of a maturity structured tumor cord cell population. Discr. Cont. Dyn. Syst. B. 2004; 4: 115-134.
5. Ambrosi D., Preziosi L. On the closure of mass balance models for tumour growth. Math Mod Meth Appl Sci. 2002; 12: 737-754.
6. Sciume G, Shelton S, Gray W, Miller C, Hussain F, Ferrari M, Decuzzi P, Schrefler B. A multiphase model for three-dimensional tumor growth. New J Phys. 2013; 15: 015005-5010.
7. Kepler TB, Elston C. Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations. Biophys J. 2001; 81: 3116-3136.
8. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977; 81: 2340-2361.
9. Bratsun D, Volfson D, Hasty J, Tsimring L. Delay-induced stochastic oscillations in gene regulation. Proc Natl Acad Sci USA. 2005; 102: 14593-14598.
10. Bratsun DA, Volfson DN, Hasty J, Tsimring LS. Non-Mark-ovian processes in gene regulation. In: Noise in Complex Systems and Stochastic Dynamics III. Proceedings of SPIE. 2005; 5845: 210-219.
11. Alarcon T. A cellular automaton model for tumour growth in inhomogeneous environment. J Theor Biol. 2003; 225: 257274.
12. Kim Y, Stolarska MA, Othmer HG. A hybrid model for tumour spheroid growth in vitro I: Theoretical development and early results. Math Mod Meth Appl Sci. 2007; 17: 17731798.
13. Bratsun DA, Merkuriev DV, Zakharov AP, Pismen LM. Multiscale modeling of tumor growth induced by circadian rhythm disruption in epithelial tissue. J Biol Phys. 2016; 42: 107-132.
14. Bratsun DA, Zakharov AP, Pismen LM. Chemo-mechanical modeling of Tumor Growth in Elastic Epithelial Tissue. In: AIP Conference Proceedings. 2016; 1760: 020007/1-11.
15. Viktorinova I, Pismen LM, Aigouy B, Dahmann C. Modelling planar polarity of epithelia: The role of signal relay in collective cell polarization. J R Soc Interface. 2011; 8: 10591063.
16. Salm M, Pismen LM. Chemical and mechanical signaling in epithelial spreading. Phys Biol. 2012; 9: 026009-026023.
17. Bratsun D, Zakharov A, Pismen L. Modeling of tumour growth Induced by circadian rhythm disruption in epithelial tissue. Emerg, Compl Comput. 2015; 14: 295-306.
18. Zakharov AP, Pismen LM. Reshaping nemato-elastic sheets. Eur Phys J E Soft Matter. 2015; 38: 75-81.
19. Zakharov A, Bratsun D. Synchronization of circadian rhythms at scale of gene, cell and whole organism. Emerg, Compl Comput. 2014; 8: 345-355.
20. Greene MW. Circadian rhythms and tumor growth. Cancer Lett. 2012; 318: 115-123.
21. Sahar S, Sassone-Corsi P. Regulation of metabolism: the circadian clock dictates the time. Trends Endocrinol. Metab. 2012; 23: 1-8.
22. Martinez-Torres C, Laperrousaz B, Berguiga L, Boyer-Provera E, Elezgaray J, Nicolini FE, Maguer-Satta V, Arneodo A, Argoul F. Deciphering the internal complexity of living cells with quantitative phase microscopy: a multiscale approach. J Biomed Opt. 2015; 20(9): 096005.
23. Friedl P., Gilmour D. Collective cell migration in morphogenesis, regeneration and cancer. Nat Rev Mol Cell Biol. 2009; 10: 445-457.
24. Pittendrigh CS. Temporal organization: reflections of a Darwinian clock-watcher. Annu Rev Physiol. 1993; 55: 16-54.
25. Lakin-Thomas PL, Brody S. Circadian rhythms in microorganisms: New complexities. Annu Rev Microbiol. 2004; 58: 489-519.
26. Gonze D, Bernard S, Waltermann C, Kramer A, Herzel H. Spontaneous synchronization of coupled circadian oscillators. Biophys J. 2005; 89: 120-129.
27. Tsimring LS. Noise in biology. Rep Prog Phys. 2014; 77(2): 026601.
28. Szekely T, Burrage K. Stochastic simulation in systems biology. Comput Struct Biotechnol J. 2014; 12: 14-25.
29. Lemerle C, Di Ventura B, Serrano L. Space as the final frontier in stochastic simulations of biological systems. FEBS Letters. 2005; 579: 1789-1794.
30. Yu EA, Weaver DR. Disrupting the circadian clock: Gene-specific effects on aging, cancer, and other phenotypes. Aging. 2011; 3: 479-493.
31. Yeoma M, Pendergasta JS, Ohmiyab Y, Yamazaki S. Circa-dian-independent cell mitosis in immortalized fibroblasts. Proc Natl Acad Sci USA. 2010; 107: 9665-9670.
13. Bratsun DA, Merkuriev DV, Zakharov AP, Pismen LM Multi-scale modeling of tumor growth induced by circadian rhythm
disruption in epithelial tissue. J Biol Phys. 2016; 42: 107132.
14. Bratsun DA, Zakharov AP, Pismen LM. Chemo-mechanical modeling of Tumor Growth in Elastic Epithelial Tissue. In: AIP Conference Proceedings. 2016; 1760: 020007/1-11.
15. Viktorinova I, Pismen LM, Aigouy B, Dahmann C. Modelling planar polarity of epithelia: the role of signal relay in collective cell polarization. J R Soc Interface. 2011; 8: 1059-1063.
16. Salm M, Pismen LM. Chemical and mechanical signaling in epithelial spreading. Phys Biol. 2012; 9: 026009-026023.
17. Bratsun D, Zakharov A, Pismen L. Modeling of tumour growth Induced by circadian rhythm disruption in epithelial tissue. Emerg Compl Comput. 2015; 14: 295-306.
18. Zakharov AP, Pismen LM. Reshaping nemato-elastic sheets. Eur Phys J E Soft Matter. 2015; 38: 75-81.
19. Zakharov A, Bratsun D. Synchronization of circadian Rhythms at scale of gene, cell and whole organism. Emerg Compl Computat. 2014; 8: 345-355.
20. Greene MW. Circadian rhythms and tumor growth. Cancer Lett. 2012; 318: 115-123.
21. Sahar S, Sassone-Corsi P. Regulation of metabolism: the circadian clock dictates the time.Trends Endocrinol Metab. 2012; 23: 1-8.
22. Martinez-Torres C, Laperrousaz B, Berguiga B, Boyer-Provera E, Elezgaray J, Nicolini FE, Maguer-Satta V, Arneodo A, Argoul F. Deciphering the internal complexity of living cells
with quantitative phase microscopy: A multiscale approach. J Biomed Opt. 2015; 20(9): 096005.
23. Friedl P, Gilmour D. Collective cell migration in morphogenesis, regeneration and cancer. Nat Rev Mol Cell Biol. 2009; 10: 445-457.
24. Pittendrigh CS. Temporal organization: reflections of a Darwinian clock-watcher. Annu Rev Physiol. 1993; 55: 16-54.
25. Lakin-Thomas PL, Brody S. Circadian rhythms in microorganisms: New complexities. Annu Rev Microbiol. 2004; 58: 489-519.
26. Gonze D, Bernard S, Waltermann C, Kramer A, Herzel H. Spontaneous synchronization of coupled circadian oscillators. Biophys J. 2005; 89: 120-129.
27. Tsimring LS. Noise in biology. Rep Prog Phys. 2014; 77(2): 026601.
28. Szekely T, Burrage K. Stochastic simulation in systems biology. Comput Struct Biotechnol J. 2014; 12: 14-25.
29. Lemerle C, Di Ventura B, Serrano L. Space as the final frontier in stochastic simulations of biological systems. FEBS Letters. 2005; 579: 1789-1794.
30. Yu EA, Weaver DR. Disrupting the circadian clock: Gene-specific effects on aging, cancer, and other phenotypes. Aging. 2011;3: 479-493.
31. Yeoma M, Pendergasta JS, Ohmiyab Y, Yamazaki S. Circa-dian-independent cell mitosis in immortalized fibroblasts. Proc Natl Acad Sci USA. 2010; 107: 9665-9670.
