PROTEIN PATTERN FORMATION INDUCED BY THE JOINT EFFECT OF NOISE AND DELAY IN A MULTI-CELLULAR SYSTEM

. We explore the combined eﬀect of the intrinsic noise and time delay on the spatial pattern formation within the framework of a multi-scale mobile lattice model mimicking two-dimensional epithelium tissues. Every cell is represented by an elastic polygon changing its form and size under pressure from the surrounding cells. The model includes the procedure of minimization of the potential energy of tissue. The protein ﬂuctuations in the tissue are driven by transcription/translation processes in epithelial cells exchanging chemical and mechanical signals. Network architecture includes a simple autorepressor model with time-delayed negative feedback, in which the only gene deﬁnes the oscillatory activity. Simultaneously, the expressed protein of the autorepressor acts as a positive regulator of the signaling protein by activating its transcription. The signaling species is assumed to spread from one cell to the other by the diﬀusion mechanism. We provide both deterministic and stochastic descriptions. The numerical simulation of spatially-extended stochastic oscillations is performed using a generalized Gillespie algorithm. We developed this method earlier to account for the non-Markovian properties of random biochemical events with delay. Finally, we demonstrate that time delay, intrinsic noise, and spatial signaling can cause a system to develop the protein pattern even when its deterministic counterpart exhibits no pattern formation.


Introduction
Since the second half of the 20th century, molecular biology has been increasingly transformed from a descriptive science to a mathematical one.The main reason why this process took so long is that the experimental technologies used in the genetic analysis were insufficient to make accurate measurements of protein fields in space and time.The isolation of the green fluorescent protein GFP, whose gene is now widely used as a tool for studying gene expression [43], has become a crucial step along this path.The ability to register, measure, and, consequently, control the protein field gave rise to a new research area of synthetic biology.The work [16] was a pioneer in synthetic biology.The authors theoretically modeled and then experimentally constructed a plasmid consisting of three different genes, which do not occur naturally in such a combination in one genome.The invented scheme generated oscillatory gene expression, which was measured experimentally using a fluorescent protein.This experiment with the repressilator proved that this could be the way to design synthetic gene D. BRATSUN networks with a prescribed dynamical behavior.Then the different effects in the synthetic gene circuits were experimentally and theoretically studied.The papers [13,18] showed that the gene circuits could do some simple arithmetic calculations.Logic operations made by the synthetic network were demonstrated in [4].Many works focused on studying the oscillations.The scholars proposed and experimented with more complicated oscillatory schemes than the repressilator's scheme [19,49].As shown in [42], the accuracy of the molecular clock of the repressilator could be improved by the simple scheme upgrade.The ensembles of interacting repressilators were studied in [52].Spatial protein waves were detected experimentally in the bacterial ensembles in [14].
It is worthy to notice that the small number of reactant molecules involved in gene regulation can lead to significant fluctuations in mRNA and protein concentrations.There have been numerous studies devoted to the influence of such noise at the regulatory level since pioneering works in the early 2000s [28,30,44].For a good review of developments in this area, see [51].
In recent years, the interest of researchers has been attracted by the time delay of gene regulation processes, as well as more general biological systems with memory [6,14,15,27,41,50,51,55].The transcription-translation processes are compound multi-stage reactions involving the sequential assembly of long molecules.Thus, these processes are distributed in space, extended in time, and proceed with some delays.Until delays are small compared with other time scales characterizing the genetic system, one can safely ignore them in simulations.However, if the delays become longer than other processes, the system has to be considered as non-Markovian, and one should account for it in both deterministic and stochastic descriptions [5,6].It is interesting to note that nature successfully uses the time delay mechanism to excite robust and accurate oscillations in terms of the period and amplitude.For example, the fact that delayed-induced stochastic oscillations can occur during transcriptional regulation was established in studies devoted to circadian rhythms in Neurospora, Drosophila, and others [15].It is widely accepted now that these oscillations are caused by delays in certain elements of gene regulation networks [46].It is plausible that the role of time delays in circadian rhythms has come to light because the delays in the corresponding reactions are particularly long (several hours) compared with other characteristic times of the system.We used circadian rhythms models to describe the carcinoma occurrence due to a local failure to synchronize rhythms in cells [9].
We studied the combined effect of the intrinsic noise and time delay on the temporal behavior during gene regulation in [5,6].We have suggested a modification of the well-known Gillespie algorithm [20] widely used to simulate statistically correct trajectories of the state of a chemical reaction network that accounts for delay [5].Based on this technique, we showed that quasi-regular fluctuations could arise in the stochastic system with delay even when its deterministic counterpart exhibits no oscillations [6].Since that papers, there have been a lot of works developing this line of research [12,27,37,41,53,55], to mention but a few.See also [25,40,48,51] for different reviews.They mainly focus on further improving the algorithm and studying the temporal dynamics of the various gene systems with delays.
Several years ago, Lemerle et al. [34] have noticed that "space is the final frontier in stochastic simulations of biological systems".The problem is that a massive amount of spatio-temporal experimental data has been accumulating, but stochastic models of biochemical processes mainly focus on temporal dynamics.If spatial stochastic simulations of Markovian processes have made considerable progress in the past years [11,34,35,38,51], the examples of studies of non-Markovian stochastic systems are still rare.For instance, Marquez-Lago et al. [36] discussed how the molecular wandering in space could be taken into account in purely temporal models through distributed delays.Danino et al. [14] have given the remarkable experimental data with the pattern formation of delay-induced rhythms in a population of E. Coli, but did not provide the stochastic modeling.
In this paper, we study the spatial effects caused by the combined action of intrinsic noise and time delay, which are simultaneously present in a multi-cellular system.To construct the latter, we use a simplified version of the chemo-mechanical model of the epithelial tissue, which first was developed in [45], then was applied to simulate the carcinoma growth [9,10] and some other problems [31].In recent years, hybrid models have become more and more popular, see, for example, the papers [2,47].Using the lattice approach allows us to avoid the development of a reliable Gillespie-like algorithm for spatial non-Markovian processes and to apply the algorithm suggested in [5,6].We use a single-gene auto-repressor model with a constant delay at a single-cell level complemented by the dynamics of the signaling species penetrating cell membranes.

Single-cell dynamics
Let us first consider the processes that take place in one cell.In this work, we assume that the dynamics of each cell are determined by a single-gene protein synthesis with negative auto-regulation (Fig. 1(a)).It is a popular motif in genetic regulatory circuits, and its temporal dynamics were analyzed within both deterministic and stochastic framework [5,6,30].This model is relatively simple but still maintains a high degree of biological relevance.
As is known, the classical Goodwin model is a minimal oscillator based on a negative feedback loop [22,23].It includes a three-step chain of reactions for the clock gene mRNA, the clock protein, and the transcriptional repressor that inhibits the production of the first component.The problem is that no limit-cycle oscillations can occur in this model until the degree of a Hill term becomes greater than 8 [24].However, a generalization of the Goodwin model to the case of time-delayed feedback is much more practical because one can obtain the oscillatory behavior for the smaller value of a Hill coefficient [21].Moreover, as noted above, a time delay is the common cause of oscillations in genetic systems due to gene regulation.The reason is that these processes are typically slow and comprise multistage biochemical reactions engaging the sequential assembly of long molecules and therefore are likely to generate time delays.Thus, we generalized and studied in [5,6] the model of a singlegene autorepressor, taking into account that the processes of transcription/translation of the protein take a finite time τ to complete.
Let us assume that the only gene x is active in the cell (Fig. 1(a)).Suppose that its protein can exist both in the form of folded monomers X and dimers X D .Both forms coexist and interact with each other via the reactions of dimerization: where k +d , k −d stand for reaction rates.
Let us denote the unoccupied and occupied state of the promoter site of the gene as D 0 and D 1 , respectively.The transitions between operator states for each protein occur with rates k +1 , k −1 when some dimer binds to the promoter or unbinds from it respectively: A decisive role in this model comes from the delay τ in the synthesis reaction: where A stands for the rate of a time-delayed production of protein monomer.One can see that the chemical state of the operator sites D∈ {D 0 ,D 1 } determines the production of corresponding protein at time t + τ .If the operator-site at time t is unoccupied (D 0 ), the protein may appear at time t + τ .Otherwise, if the operator is occupied (D 1 ), the production at time t + τ is blocked.Thus, we have set a negative feedback loop in this model.The more the system synthesizes the protein, the more dimers appear, and the more likely one of the dimers will bind to the operator site and inhibit transcription.In the model, we ignore the appearance of mRNA as an intermediate form in synthesizing a protein since we believe that the time delay formally accounts for the translation process.Finally, we take into account the effect of protein degradation with a rate of B: To describe the effect of intercellular signaling, we introduce the signaling species T positively regulated by the X monomers.For simplicity, let us assume that the signaling protein T is synthesized without delay, and its expression is regulated by X monomers, the simplest form of the autorepressor protein.Let the chemical state of the operator site D T ∈ {D T 0 ,D T 1 } of the tr-gene determines the production of the T protein (Fig. 1(a)).If the operator is occupied (D T 1 ), the protein may be produced immediately with a certain probability A T in a unit time.Otherwise, if the operator is unoccupied (D T 0 ), the production of signaling protein T is blocked.The transitions between operator states occur with rates k +2 , k −2 , are Taking into account the positive feedback between X monomers and the expression of the signaling protein, we write the synthesis reaction in the following form: where A T is the reaction rate.
Finally, we assume also that once a signal has come in some cell, it is converted back into X protein with the rate B T .The activation is assumed to be linear in the signal concentration: (2.7) Thus, the reactions (2.6), (2.7) have the positive feedback with (2.5) through the reaction rate in (2.6).The protein T does not affect the processes of transcription/translation of x-gene.Altogether, equations (2.1)-(2.7)define the kinetics of gene regulation at a single-cell level.

Deterministic description
Within the framework of deterministic analysis, we neglect any kind of noise, considering the influence of both intrinsic and extrinsic noise on the dynamics of the system to be insignificant.Following the standard procedure, we obtain the following set of delay-differential equations (DDEs) derived from the reaction equations (2.1)-(2.7): where the subscripts refer to cells.Here, the average concentrations of the X monomers, the X D dimers and the T monomers in i-th cell with the area S i (t) at time t are denoted as x(t), y(t), and θ(t), respectively.We should notice that our cellular model, as will be shown below, is dynamic.So, the concentration of species can change not only due to chemical kinetics (2.1)-(2.7)but also due to changes in the area of a cell under external pressure in the tissue.In (3.1)-(3.7), the continuous variables s 0i (t), s 1i (t), s T 0i (t) and s T 1i (t) stand for the average number of unoccupied and occupied operator sites of x and tr genes in i-th cell at time t respectively.The last term in equation (3.7) means that chemical signals between cells are delivered via a diffusion mechanism.We assume that the T protein is transported diffusively from one cell to the other.This flux does not depend on the distance between the two cells i and j but is proportional to the boundary length L ij and the concentration difference θ j (t) − θ i (t).The intensity of the process is also determined by the diffusion coefficient α.Generally, it implies that the intercellular transport includes only the transfer through cell membranes.The link between sub-cellular and macroscopic scales is established through the equation (3.7) since the θ field is global for the whole tissue.
In addition to equations (3.1)-(3.7),we assume that where the integer constant on the right side of each equation stands for the copy number.It means how many copies of the monomer can be synthesized in a single act of transcription/ translation.For simplicity, we assume that the copy numbers in both cases are equal to 1.
The main approximation we make here is an assumption that the dimerization reactions (2.1) and binding/unbinding (2.2), (2.5) are fast in comparison with production/degradation reactions (2.3)-(2.7): (3.9) The assumption (3.9) makes us divide the dynamics of the system into slow and fast segments.The final purpose of the analysis is to arrive at the motion equations on the slow manifold in phase space.Generally, a sequential averaging procedure can be done by multiple-scale analysis, but it can be simplified for the set of DDEs (3.1)-(3.7)by giving a zero value to the respective derivatives in the left parts of the equations.So, we suppose that dynamics of operator-sites and dimers quickly enter into a local equilibrium, in which the concentrations of species satisfy where By solving (3.8) and (3.10) together for each operator-site, we arrive at the expression for average values of open operator-sites: To close the equation system, one should write the equation for the total number of X protein molecules x all in a monomer form, where the first term in the right part determines the number of free X monomers; the second term defines the total number of free X dimers (two monomers in one dimer molecule); the third term is responsible for the number of dimers bound with the operator-site of the x gene; finally, the last term gives the number of dimers bound with the operator-site of the tr gene.In general, the third and the fourth terms are small compared to the first and the second terms, and, therefore, we can neglect them.The final equation of the autorepressor dynamics is obtained by differentiating equation (3.12) with respect to time: Here, we neglect the second-and higher-order expansion terms.By taking into account (3.10)-(3.12), the system of DDEs (3.1)-(3.7)can be reduced to the following form: The prefactor appearing in equation (3.13) does not affect the possible equilibrium states of the system but significantly changes their stability and nonlinear dynamics.At a single-cell level without the intercellular signaling (θ i = 0, A T = 0), one can prove that equations (3.13), (3.14) have only a positive steady-state solution.In the case of ε 1, δ 1 one can derive the following approximate solution: The assumption that the constants are small ε 1, δ 1 is quite reasonable since, for the stability of the process, the rate of the reverse reaction must be an order of magnitude higher.As a result, equilibrium is determined by the exact balance of protein production and degradation: x * ≈ A/B.One can show that the fixed point (3.15) in the absence of delay is absolutely stable, which reflects the precise work of the autorepressor: negative feedback blocks the operator site just enough to maintain the desired level of gene expression.The prefactor, which appears due to the correct counting of protein molecules in the system (3.12),does not affect the equilibrium states but significantly affects the stability of the steady-state and the phase dynamics.
By linearizing equations (3.13)-(3.14)near (3.15), and looking for a solution of the form x(t) ∼ e λt , where λ = χ + iω, we obtain closed-form analytical solutions for λ: ) where W (z) stands for the Lambert function defined as W (z)e W (z) = z.First of all, let us suppose that λ = 0 in (3.16), (3.17).In this case, we can get the following expression: where we used equation (3.15).As far as all left values in (3.18) are positive, this expression can not be satisfied as it is.It means that there is no bifurcation in simple eigenvalue for steady-state x * , and it is absolutely stable with respect to monotonic disturbances.The neutral curve for the Andronov-Hopf bifurcation is determined by the condition λ = Iω (Fig. 2(a)).Here, ω is the oscillation frequency of the limit cycle at the moment of its branching.Near the bifuraction point, it approximately equals to ω * ≈ π/τ .Figure 2(b) shows the power spectrum of the periodic oscillations calculated in the instability domain (it is indicated by the upper black point in Fig. 2(a)).
Since there is no positive feedback loop in a single-cell model (2.1)-(2.7),and the internal source of energy of the system due to the processes of gene transcription/translation is balanced by regulation through negative feedback, the system does not demonstrate anything more complicated than the limiting cycle in the supercritical region of parameters.The oscillation arises due to the delay and reflects the inability of the autorepressor to respond in time to changes in its state.It is the common reason for the development of oscillatory behavior in the control theory.

Stochastic description
To describe the spatial stochastic effects, we use a hybrid model constructed as follows.The dynamics of the X protein in each cell is obtained by performing direct Gillespie simulations of a single-gene auto-repressor model given by the reactions (2.1)-(2.7).The modified version of the Gillespie algorithm was developed in [6].This algorithm accounts for the non-Markovian properties of random biochemical events with delay.Generally, Earlier, we demonstrated in [5,6] that the time delay coupling with intrinsic noise could cause a genetic system to be oscillatory even when its deterministic counterpart exhibits no oscillations.But what about space?
In stochastic simulations we assume values X i , X Di , T i , D 0i , D 1i , D T 0i , D T 1i , to be integer.The algorithm then contains the following steps for each i-th cell: 1. Input values for initial state of channels X i , X Di , T i , D 0i , D 1i , D T 0i , D T 1i , set t = 0. 2. Compute propensities a µi , µ = 1..7. 3. Generate two uniform random numbers u 1i , u 2i .4. Compute the time step until the next reaction: t i stoch = −(ln u 1i )/ µ a µi . 5. Check if there has been either a delayed reaction scheduled at time t delay or a diffusion event scheduled at time t dif f to occur first within the range between t and t + t i stoch : A) If there has been scheduled a delayed reaction at time t delay < t dif f to occur first in the range between t and t + t i stoch : a) If yes, then the steps 2-4 are ignored, time advances to t = t delay , the time-delayed species is updated.Proceed to step 2. b) If no, go to the step B. B) If there has been a diffusion event should occur at time t dif f within the range between t and t + t i stoch : a) If yes, then the T i species is updated according to following finite-difference formula: where [...] stands for the integer part of the expression, S i is the area of i-th cell.
Update the positions of the lattice nodes and mechanics of the system.Proceed to step 6. b) If no, go to the step 6. 6. Find the channel of the next reaction: take µ to be the integer for which , where a 0i is the total propensity, a 0i = M ν=1 a µi .7. Update time t → t + t i stoch .If the selected reaction µ is common, update the species.If the reaction is delayed, the update is postponed for time t + τ .
Proceed to step 2.

Chemo-mechanical model of two-dimensional tissue
Protein fields arise in the tissue of living organisms.The tissue consists of a huge number of cells exchanging both chemical and mechanical signals.Cells differ in shape, size, physical properties, and phenotype.For example, cells of the mesenchymal phenotype can move through the tissue since they are not attached to other cells by desmosomes.On the other hand, cells of the epithelial phenotype can divide.All of the above can generate mechanical waves in the tissue controlled by excess pressure.Generally, tissue cells try to take such a shape, in which the total potential energy would be less.Therefore, it seems almost obvious that the simulation of the protein dynamics should be carried out not on a uniform static grid but in a movable lattice that mimics the behavior of living tissue.
Epithelium can be defined as a relatively avascular aggregation of cells, which are in apposition over a large part of their surfaces, and are specialized for absorptive, secretory, protective, or sensory activities.Epithelia may occur as sheets of cells, characterized as a covering and lining epithelium, or as solid aggregations of cells, as in glandular organs.The chemo-mechanical model we use focuses on the first type and presents the epithelium as an elastic two-dimensional array of cells, approximated by polygons, and includes the simulation of the dynamics of separate cells adapting to environmental stresses.The model has been developed initially in the work [45], and then has been applied to simulate the carcinoma cancer in [9,10,31].
The evolution of the system starts from hexagonal cells constituting the plane.In the course of proliferation and division, the initial lattice deforms, incorporating irregular polygons.But, generally, the system was calibrated in such a way, as the hexagonal cell is the most probable form of a cell.Also, we neglect the curvature of the layer (presumed small compared to the cell size) and thickness inhomogeneities.
We use a simplified version of the model with the following set of properties: -allowing for changing the size and shape of cells; -allowing for tissue spreading via the mechanism of cell division (Fig. 3  The mechanical model is based on the elastic potential energy U of the tissue, defined by summing up the contributions of the perimeter L and the area S of each cell [9,17,45]: Here the coefficient κ is attributed to the action of active contractile forces within the cell cortex, η is the elastic constant that reflects resistance to stretching or compressing the cell vertically when its area decreases or increases with respect to the reference cell area S 0 .Vertices of the polygons representing the cells form a mobile lattice.The evolution of the tissue proceeds by moving the lattice nodes indicated by blue circles in Figure 3.We define the mechanical force acting on any jth node in the following way: where R j denote the position of the node and F st j is an uncorrelated zero-mean stochastic input slightly disturbing the system.
We consider the internal movement of cells in the tissue as a strongly overdamped process, so the equation of motion can be written as follows: where V i is the velocity of ith node, K is the mobility coefficient.The threshold force F 0 has been introduced in (3) to take into account the situation when the node remains immobile even if the force is not zero (see the discussion in [45]).H stands for the Heaviside step function.
An important feature of the model is the ability of cells to divide.It 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. 3(a)).To minimize the growing disorder in the distribution of nodes in the process of cell division, the probability to divide p was connected to the number of vertices of the polygon n j according to the following formula (see [10,45] for more detail): where q is a distribution constant and p 0 is a scaling factor.One can see that for q > 1, the polygons with more than 6 sides are more often divided.Another mechanism that increases the liquidity of the tissue and excludes the severe deformation of cells is intercalation [10,26,45].We introduce here a special parameter l 0 that defines the cases in which intercalation can occur.If the length of the border separating two cells becomes less than l 0 , it is substituted by a link of a slightly larger length in the normal direction.Intercalation is known to be important in many tissue reshaping processes.Altogether, equations (5.1)-(5.4)define the mechanical dynamics of tissue on the cellular level.
As is known, living tissue is a typical chemo-mechanical elastomer.It means that the chemical signals generated by gene regulatory circuits in cells directly affect tissue elasticity.Conversely, mechanical changes in tissue can affect the biochemical processes of gene regulation in cells.In this work, we consider a simplified model, in which chemistry does not affect cell properties.However, the reverse effect of tissue mechanics on gene regulation is maintained through changes in the shape and area of the cell (Fig. 3(c)).
In our model, a cell is a pixel of spatial resolution.We ignore any heterogeneity that occurs within a cell.However, when changing the cell area, we should take into account that the species concentrations in the cell also change, even while the number of molecules remains the same.Figure 3(c) illustrates how the flux of the signaling protein arises due to a sharp decrease in the size of the left cell.Equation (4.1) takes this effect into account when calculating protein fluxes.Thus, our algorithm includes the calculation of signal protein fluxes between cells based on the calculation of concentrations taking into account changes in the shape and size of cells.

Numerical results
To study the spatial effects numerically within the deterministic description, the set of DDEs (3.13), (3.14) was solved using the explicit Euler method.Its stability was warranted by a sufficiently small time step ∆t chem = 0.005.The time step for the calculation of the molecular processes in cells was synchronized with the step ∆t mech = 0.01 of simulation of the mechanical evolution of the tissue.To avoid storing a large amount of data for long time delay intervals, we have used a numerical simulation algorithm [7] based on adaptive storing in the computer memory of some selected nodal data and a subsequent interpolation to calculate the intermediate values.This adaptive numerical scheme for data storing can speed up the calculation up to 8 times without visible deviations from benchmark simulation.In the case of stochastic simulations, we used the hybrid algorithm described above.
The initial configuration of the system is a hexagonal lattice comprising 1560 cells filing the the rectangular domain −30 ξ 30, −35 η 35.The typical values of the parameters governing the tissue mechanics hereinafter are as follows: κ = 1.0, η = 1.0, S 0 = 3 √ 3/2, F 0 = 0.02.In the case of deterministic description given by equations (3.13), (3.14), the numerical simulation in each cell start from a random value of the initial concentration x(0).In the case of stochastic simulation, at the very beginning, we set a random number of X monomers, an open operator site D 0 = 1, and a zero number of dimers Y = 0. Then the initial phase of oscillations is determined by the X protein level.As is known from experimental observations of the processes of gene regulation, the dynamics of the operatorsite and the processes of dimerization proceed an order of magnitude faster than the processes of protein synthesis and degradation (3.9).If the value of the characteristic time for the former processes is about seconds-minutes, then for the latter -hours.We used the assumption (3.9) when deriving a closed-form model (3.13), (3.14) providing the deterministic description of the system slow dynamics.Therefore, in order to make a correct comparison of stochastic simulation results with those of the equations (3.13), (3.14), we must not only set the ratio of the rates of fast reactions ε and δ, but also control in simulations the reaction rates k i themselves.In this work, we used the following values k +d = 200, k −d = 1000, k +1 = 100, k +2 = 100, k −1 = 1000, and k −2 = 1000, which satisfy (3.9).
Figure 2(c, d) present the power spectra of a time-delayed stochastic signal generated on a single-cell level respectively below and above the neutral curve.These points on the stability map are indicated in Figure 2(a) by the black squares.One can notice a remarkable row of peaks indicating that some periodicity in the signal should exist.The first peak corresponds to the frequency ω * = π/τ ≈ 0.524.Also, one can observe the strong harmonics of the fundamental frequency given by ω n = (2n − 1)ω * .Thus, we can conclude that the time delay coupled with the intrinsic noise can force a system to be oscillatory even when its deterministic description predicts no oscillations.
Before we proceed to the description of the numerical results for a spatially extended system, let us briefly figure out how the system evolves.One must keep in mind that there are two evolving subsystems.The first one includes the mechanical evolution of the epithelium governed by deterministic equations (5.1)- (5.4).The main cause of changes in cell tissue is cell division, which is given by a random probability distribution (5.4).Cell division changes the local potential energy (5.1) and forces the cells to slightly change in shape, size, and position in space.The second subsystem includes the processes of gene regulation in cells.The deterministic dynamics is   2(a)).The frames from left to right and from up to down correspond to times 260, 296, 332, 368, respectively.These consecutive times are separated by 6τ to produce the stroboscopic effect, which is represented by a instantaneous samples at a sampling rate close to 3 periods of the oscillations.The frequency of base background oscillations is ω * ≈ π/τ (the period is about 2τ ).The initial configuration of the cellular tissue is a regular hexagonal lattice comprising 1560 cells.The domain of the integration is −30 ξ 30, −35 η 35.The colorbar indicates the mapping of protein concentration values into the colormap.described by equations (3.13), (3.14).Numerical simulation of stochastic dynamics implies finding the number of protein molecules in each cell separately following the algorithm described in Section 4. The calculations of intercellular diffusion and the changes in the epithelial tissue are done on the same uniform time grid.At each time step, we calculate the signal protein fluxes that provide cell-to-cell communications (see Eq. (4.1)).We have to take into account that the time step of stochastic simulations in each cell changes randomly following the algorithm.To calculate the intercellular flux of signaling protein, we use the concentration values established in a given cell at the last step of stochastic dynamics.
Figure 4 presents the results of numerical simulation of a spatially-extended system made in the framework of the deterministic description.The parameter values are taken above the Hopf bifurcation curve (this point is indicated in Fig. 2(a) by a black square).One can see that each cell demonstrates oscillatory behavior.In the figure, the X protein field is shown for six consecutive times taken within one period of oscillations.The nonlinear dynamics include the slow development of spiral traveling wave pattern, which arises against the synchronized field oscillating in the background.The waves develop slowly due to the negative feedback of the autorepressor.
Let us consider two examples of spatial stochastic simulations.Figures 5 and 6 presents the time evolution of the stochastic pattern formed by the X monomers distributed in the tissue (see also GIF-animations in Supplementary materials).The values of governing parameters are the same as in Figure 4. We found that the nonlinear dynamics of a spatially-extended system consists of two distinct oscillatory modes, just like it was in the deterministic case.One is a quasi-standing wave pattern oscillating with the frequency ω * (see Fig. 5(a)).The second oscillatory mode consists of traveling waves, which arise from selected initial disturbances (see Fig. 6).In fact, the stochastic pattern looks very similar to its deterministic counterpart obtained for the same parameters.Small differences in the pattern formation (for example, the wavelength of the structure in Fig. 6 looks a little larger) within the deterministic and stochastic descriptions can be explained by the incomplete equivalence of these approaches.We should keep in mind the set of DDEs (3.13), (3.14) was obtained under condition (3.9). Figure 5(b) shows the temporal evolution of the distribution of tissue cells over the oscillation phase, which is defined by the number of the X monomers in the cell.Each frame in Figure 5(b) corresponds to a frame in Figure 5(a).We can notice that the system evolution within one oscillation period includes the stage of activity and passivity of the protein field.In the first case, the most probable phase in the distribution corresponds to a high protein level in cells, and in the second, the most probable protein level drops sharply.More importantly, at the moment when the system demonstrates clear clustering of cells by oscillation phases, the distribution acquires two maxima (see, for example, the first frame in Fig. 5(a)).Thus, we can conclude that a two-humped distribution structure is a marker of the clustering effect in tissue.It is worthy to note that the clusters demonstrate oscillatory behavior.At certain moments, the cell phases align since the protein level rises in one part of the cells and falls in the other, and the two-humped distribution structure disappears.
As is known, the clustering in the system with a large number of elements exchanging chemical signals has become at the center of attention of many scientists recently (see, for example, [32]).It is widely believed now that clustering results in the differentiation of cells in organs.
Let us consider further the case of parameter values far below the Andronov-Hopf bifurcation curve indicated in Fig. 2a).The system starts with a random distribution of X protein in cells and quickly falls into a fully synchronized mode of oscillations with the frequency ω * .Besides, we found the clustering effect when the cells form two approximately equal communities, which oscillate in anti-phase.In Figure 7(a), this effect clearly manifests after a sufficiently long integration.Figure 7(b) shows the dynamics of the distribution of cells by the concentration of X protein.To enhance the effect, we chose one frame of evolution at time t = 422 (see Fig. 7), split the integration region into two equal parts (ξ < 0 and ξ > 0), and separately calculated the cell distribution for the two subdomains.Figure 8 presents two histograms for subdomains.It is now clearly seen that there are two clusters of cells in the tissue.Finally, Figure 9 shows the dynamics of the cell clusters over time.
We also found that the clustering effect does not occur for any values of the control parameters in the region below the neutral curve.The further the parameter values are from the oscillation excitation curve, the longer it takes to obtain clustering.The effect is enhanced by increasing the positive feedback between the X protein and the signaling protein.One can do it, for example, by setting the copy number in (3.8) to be greater than 1.
Finally, we found that small changes in the structure of cellular tissue associated with intercalation and division do not fundamentally affect the properties of the protein field.However, in the case of a more drastic change in the tissue mechanics, this may not be the case.For example, a sudden local tissue compression may lead to a sharp increase in the protein concentration in cells being in the compression zone.We did not consider these effects in this work.

Discussion
In our opinion, the main lesson that we can learn from this study is that the processes of gene regulation occurring in tissues in vivo never stop and proceed unsteadily.One can develop a mathematical model, which predicts the absolute stability of the system for some values of control parameters.However, the real genetic system is always noisy.The intensity of intrinsic noise increases sharply due to the small number of reactant molecules involved in gene regulation.There are also extrinsic sources of noise, which are attributable to a noisy cellular environment.All this together leads to significant fluctuations in intracellular mRNA and protein concentrations.Also, we should keep in mind that the transcriptional and translational processes are multistage   7(a)).The blue and red histograms represent the statistics for cells from the subdomains ξ < 0 and ξ > 0, respectively.reactions involving the sequential assembly of long molecules.It implies that these processes must occur with some characteristic delay time.In this work, based on the autorepressor model, we demonstrated that the combined effects of time delay and intrinsic noise on gene regulation produce oscillations and weak patterning in space even when the deterministic description of the system shows only stationary behavior.
Generally, the autorepressor without delay is a common genetic control motif consisting of a gene that encodes a protein repressing its expression.This architecture has been well-studied in various contexts, and is employed ubiquitously to control the expression of repressors in both prokaryotes and eukaryotes.Several variants of the autorepressor have been experimentally realized (see, for example, papers [1,3,29] for details).
Let us discuss possible ways of experimental observation of the effect considered in the present work.The study of the effects of delay in the processes of transcription/translation was stimulated primarily by the experimental discovery in some organisms demonstrating long delays in the expression of genes responsible for the circadian rhythm in cells [8,33,46].However, we should note that the circadian rhythm problem is similar but not a direct prototype of the model we consider in this paper.First, the network architecture responsible for the circadian rhythm is more complex, as it contains positive feedback.Secondly, the synchronization of circadian rhythms in the cells of higher organisms does not occur locally due to the suprachiasmatic nuclei located in the hypothalamus.As a consequence, any pattern formation cannot occur in such a medium [21].If the external control disappears for some reasons [39,54], the weak coupling of local oscillators may result in the possible formation of collective patterns.
Instead of looking for a natural prototype of the model under the consideration, an experimenter could artificially design the gene circuit.As is known, synthetic biology has developed a new research paradigm that includes the theoretical design and experimental implementation of artificial gene circuits with prescribed dynamic behavior.These circuits can then be incorporated into living cells using vector technologies, and their dynamic behavior can be registered using fluorescent proteins.Thus, a researcher got the opportunity to design any simplified combinations of the dynamic behavior of just a few genes.The technique demonstrated in [14] is an example of how to organize experimental verification of the effect reported in this paper.
To finish this discussion, notice that we applied the simplest model, the autorepressor with time delay, demonstrating biologically relevant periodical oscillations.The prototype of such oscillatory behavior in nature can be circadian rhythms in cells.As an alternative, synthetic biology technologies could be applied to implement this scheme in practice.

Conclusion
As is known, the noise during gene expression comes about in two ways.The inherent stochasticity of transcription/translation generates intrinsic noise.The extrinsic noise refers to variation in identically-regulated quantities between different cells.In this paper, we have focused on the spatial effects of intrinsic noise.Since any living matter consists of cells, it is practical to use numerical methods based on a computational lattice mimicking a cellular compartment, such as a membrane or the interior of some part of a cell.In this paper, we have used such a lattice approach suggesting the model which accounts for mechanical interactions and chemical signal exchanges between cells.
The principal elements that we included in the model are -the minimum number of genes involved in the regulation in one cell; -the gene circuit has to demonstrate oscillatory behavior at some values of control parameters; -the delay time must be larger than all other characteristic times in the system; -the signaling or transport protein providing cell-to-cell communications; -two-dimensional epithelium-like tissue.
A new result of this study is an insight into how the excitation of quasi-regular delay-induced fluctuations of the protein concentration found in the subcritical region of the stability map [6] manifests itself in space.In this work, we used a simple single-gene autorepressor model to see what happens in space.We show that, above the Andronov-Hopf bifurcation, one observes the traveling wave pattern similar to that in the deterministic case.But the more interesting result is found below the bifurcation curve where the system within deterministic description exhibits absolute stability.As it turned out, a spatially-extended system can demonstrate not only synchronization of the protein field but also a more complicated phenomenon of clustering of cells.Thus, the combined action of intrinsic noise and feedback delays can trigger patterning.

Figure 1 .
Figure 1.(a) The network architecture of a single-gene autorepressor model with time-delayed negative feedback, in which the only gene defines the oscillatory activity.Its protein molecules indicated by red circles acts as a positive regulator of the signaling protein (blue circles) by activating its transcription.(b) The signaling species is diffused from one cell to the other within the lattice governed by the chemo-mechanical epithelium model.
(a)); -allowing for the motion of cells by the mechanism of intercalation (Fig. 3(b)); -allowing for exchanging of chemical signals between adjacent cells through their shared borders (Fig. 1(b)).

Figure 3 .
Figure 3. Elements of chemo-mechanical model of epithelium: (a) -cell division; (b) -cell intercalation; (c) -effect of cell area on cell-to-cell communications.

Figure 4 .
Figure 4. Temporal evolution of the X protein pattern in the epithelium obtained within the deterministic description given by equations (3.13), (3.14), with parameters A = 500, A T = 500, B = 5, τ = 6, α = 0.05, ε = 0.1, δ = 0.2, and σ = 0.1 taken above deterministic Hopf bifurcation curve (this point on the stability map is indicated in Fig. 2(a)).The frames from left to right and from up to down correspond to consecutive times 256, 258, 260, 262, 264, 266, respectively.These times characterize one period of base background oscillations with the frequency ω * = π/τ .The initial configuration of the cellular tissue is a regular hexagonal lattice comprising 1560 cells.The domain of the integration is −30 ξ 30, −35 η 35.The colorbar indicates the mapping of protein concentration values into the colormap.

Figure 5 .
Figure 5. (a) Time evolution of the X protein pattern in the epithelium obtained within the stochastic description with parameters A = 500, A T = 500, B = 5, τ = 6, α = 0.05, k +d = 200, k −d = 1000, k +1 = 100, k +2 = 100, k −1 = 1000, and k −2 = 1000 taken above deterministic Andronov-Hopf bifurcation curve (this point on the stability map is indicated in Fig. 2(a)).The frames from left to right and from up to down correspond to consecutive times 256, 258, 260, 262, 264, 266, respectively.These times characterize one period of base background oscillations with the frequency ω * = π/τ .The initial configuration of the cellular tissue is a regular hexagonal lattice comprising 1560 cells.The domain of the integration is −30 ξ 30, −35 η 35.The colorbar indicates the mapping of protein concentration values into the colormap; (b) Time evolution of the distribution of cells by oscillation phases defined by the number of X monomers.The frames from left to right and up to down correspond to frames shown in Figure 5(a).

Figure 6 .
Figure 6.Time evolution of the X protein pattern in the epithelium obtained within the stochastic description with parameters A = 500, A T = 500, B = 5, τ = 6, α = 0.05, k +d = 200, k −d = 1000, k +1 = 100, k +2 = 100, k −1 = 1000, and k −2 = 1000 taken above deterministic Hopf bifurcation curve (this point on the stability map is indicated in Fig. 2(a)).The frames from left to right and from up to down correspond to times 260, 296, 332, 368, respectively.These consecutive times are separated by 6τ to produce the stroboscopic effect, which is represented by a instantaneous samples at a sampling rate close to 3 periods of the oscillations.The frequency of base background oscillations is ω * ≈ π/τ (the period is about 2τ ).The initial configuration of the cellular tissue is a regular hexagonal lattice comprising 1560 cells.The domain of the integration is −30 ξ 30, −35 η 35.The colorbar indicates the mapping of protein concentration values into the colormap.

Figure 7 .
Figure 7. (a) Time evolution of the X protein pattern in the epithelium obtained within the stochastic description with parameters A = 20, A T = 20, B = 5, τ = 6, α = 0.05, k +d = 200, k −d = 1000, k +1 = 100, k +2 = 100, k −1 = 1000, and k −2 = 1000 taken below deterministic Andronov-Hopf bifurcation curve (this point on the stability map is indicated in Fig. 2(a)).The frames from left to right and from up to down correspond to consecutive times 420, 422, 424, 426, 428, 430, respectively.These times characterize one period of base background oscillations with the frequency ω * = π/τ .The initial configuration of the cellular tissue is a regular hexagonal lattice comprising 1560 cells.The domain of the integration is −30 ξ 30, −35 η 35.The colorbar indicates the mapping of protein concentration values into the colormap; (b) Time evolution of the distribution of cells by the concentrations of X protein.The frames from left to right and up to down correspond to frames shown in Figure 7(a).

Figure 8 .
Figure 8.The distribution of cells by the number of X protein in two subdomains at time t = 422 (see Fig.7(a)).The blue and red histograms represent the statistics for cells from the subdomains ξ < 0 and ξ > 0, respectively.

Figure 9 .
Figure 9.Time evolution of the distribution of cells by the concentrations of X protein obtained far below deterministic Andronov-Hopf bifurcation (this point is indicated in Fig. 2(a)).Two periods of spatially synchronized oscillations are shown.The simulation parameters are the same as in Figure 7.