{'fraction_non_alphanumeric': 0.03578805783734448, 'fraction_numerical': 0.009991833597540472, 'mean_word_length': 4.265048052604957, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Macroscopic current generated by local division and apoptosis in a minimal model of tissue dynamics

We consider a minimal computational model of tissue dynamics with two active ingredients: local particle division and apoptosis. We neglect other non-equilibrium effects such as self-propulsion. We simulated the steady state dynamics inside an asymmetric channel and we found a net macroscopic current along the channel. Although such macroscopic current in a similar geometry has been detected in swimming bacteria, our results showed that local division and apoptosis are sufficient to generate a macroscopic current, without any need for a self-propulsion/swimming mechanism. Our results might have applications in tissue engineering such as controlling tissue growth via a geometrically non-uniform substrate.


I. INTRODUCTION
It has been shown that swimming bacteria have the ability to perform useful macroscopic work. For example, when we place an asymmetric cog inside a bath full of swimming bacteria, the bacteria can rotate the cog in one direction persistently [1] (the size of the cog is around 30 times larger than the bacteria). Individually, the bacteria swim in a completely random direction. However, interactions between the bacteria and the surface of the cog can break time reversal symmetry, giving rise to a macroscopic current circling around the cog. This feature is of course not possible in an equilibrium system like passive Brownian particles. Another manifestation of the same phenomenon happens when we place swimming bacteria inside a long asymmetric pipe (similar to one depicted in Fig. 1). Experiments [2] showed that the asymmetry of the pipe induced the bacteria to move, on average, in one direction. Furthermore, it has also been shown theoretically [3] and numerically [4][5][6] that the random swimming direction (for any generic microswimmers) can be rectified through the presence of an asymmetric potential. This result is independent of any hydrodynamic or alignment interaction.
Another route to a large macroscopic current in swimming bacteria is through hydrodynamics. The rotation of the flagella around the bacterium stirs the fluid around it and a spontaneous symmetry breaking in the average orientations of the bacteria can give rise to a macroscopic flow (spontaneous flow transition [7]). This can happen in a symmetric/asymmetric potential, although the mechanism relies on hydrodynamic and alignment interactions [8,9]. Such fluids are often called active extensile or contractile fluids [10].
In this paper, we consider a minimal computational model of living tissues which can divide (mitosis) or die (apotosis). We show that local cellular division and apoptosis are sufficient to give rise to a macroscopic current inside an asymmetric channel, without any recourse to self-propulsion/swimming (see Fig. 1 and 5). Since division and apoptosis are universal properties of all living  Figure 1. The geometry considered in our model. We consider a two-dimensional box of size 2L × 2L with periodic boundary conditions in the x-direction. Along y-direction, the particles are confined in between two corrugated walls with periodicity L, peak height h and skewness ∈ [0, L]. systems, our results indicate that the ability for all living systems to perform useful macroscopic work is probably universal as well. Finally, our results might also have applications in tissue engineering, for example, using ratchet-shaped scaffold or substrate to speed up tissue growth in one direction.

II. COMPUTATIONAL MODEL
We consider a minimal computational model of living tissues in two dimensions with two key active ingredients: division and apoptosis. We shall neglect hydrodynamic interactions and any other non-equilibrium process such as self-propulsion. The aim of this paper is to demonstrate whether we can get a macroscopic current through rectification of local division and apoptosis alone.
We approximate the cells as soft circular particles which can divide into two daughter particles or die (i.e. removed from the system) [11][12][13]. Other models of tissues dynamics such as the vertex model [14] has also been considered to take into account of the confluency of biological tissues. However, recent studies [12] showed that there is no significant statistical difference between the vertex model and the particle model such as the one we use in this paper. The particles are confined inside a 2L × 2L box with periodic boundary conditions in the xdirection and asymmetric hard walls at the top and the bottom (see Fig. 1). The corrugated walls have periodicity L, peak height h and skewness value ∈ [0, L], see Fig. 1. Obviously, 2(L − h) has to be larger than the particles' diameter.
We define r i (t) = (x i (t), y i (t)) to be the centre-of-mass position of particle i, with i = 1, 2, ..., N (t), where N (t) is the total number of particles. Note that N (t) is not conserved since the particles can divide or die. We take all N (t) particles to have the same diameter d. The equation of motion for each particle i is given by the overdamped dynamics: where ζ is the friction coefficient, F ij is the force exerted on particle i by the neighbouring particle j, and F iw is the force exerted on particle i by the walls. In this model we neglect the inertia and thermal fluctuations of the particles, which are justified since the typical Reynolds number of the cells will be of order 10 −3 and cells' diameter will be of order of 10-100µm [15]. We solve Eq. (1) numerically using Euler update in time with timestep ∆t = 10 −3 . For particle-particle interactions, a soft elastic repulsive potential is used, such that the force acting on particle i due to the neighbouring particle j is given by: where k is the stiffness, r ij = |r i − r j | is the distance between particle i and particle j and H(x) is the Heaviside step function such that H(x) = 1 if x ≥ 0 and H(x) = 0 if x < 0.r ij is the unit vector pointing from r j to r i . Note that the particles only exert a repulsive force on each other if they overlap, i.e. r ij < d. Now for particle-wall interactions, the same repulsive potential is used, except that the stiffness is made much larger relative to that of particle-particle interactions. This is done to approximate the hard walls at the top and bottom of the simulation box. The resulting force acting on particle i by one side of the walls is then given by: where µ is the relative stiffness of the walls (in our simulations, we use µ = 24). r iw is the perpendicular distance For corner particles, the force from each side of the walls are added up. (c) If the particle is in contact with the edge of the wall, the direction of the force points from the edge of the wall to the centre of mass of the particle ri. This only applies if ri lies in between the two dotted lines. Figure 3. (a) Each particle i divides with rate αi, which depends on the number of contacts particle i makes with its neighbour, see Eq. (4). When the particle divides, a new particle is created at a distance ε d from the parent particle with random orientation. (b) Each particle i also dies with rate β/N . When the particle dies, it is removed from the system. from the centre of mass of the particle r i to the wall surface.r iw is the unit vector pointing from the wall to r i , in the direction perpendicular to the wall [see Fig. 2

(a)].
Again the particles only experience repulsive forces from the walls only if they overlap with the walls. In the case of a particle in contact with two sides of the walls (corner particles), the wall force from each side of the walls are then added up [see Fig. 2 If this is satisfied, the direction of the force is then computed from the edge of the wall to r i , otherwise, the rule of Fig. 2(a) applies.
We also implement particles' division and death as follows. First, each particle i divides with rate α i , which is defined to be [11]: where α 0 is the rate constant and z i is the contact number, i.e. number of contacts particle i makes with its neighbours. Overlap with another particle is counted as 1 contact while overlap with one side of the walls is counted as 2 contacts. Overlap with the wall edge [such as one in Fig. 2(c)] is counted as 1 contact. Note that the division rate goes to zero when the contact number is larger than or equal to z max . In this work, we fix z max = 6, which corresponds to the hexagonal packing/maximal packing in two-dimensions. When particle i divides, we simply add a new particle at a distance d from the parent particle in a random direction [see Fig. 3(a)]. (In our simulations, we set = 0.05d.) The daughter and the parent particles then quickly separate out in the subsequent timesteps due to strong repulsion between them. Finally, each particle i also dies with rate β/N (t). When the particle dies, we simply remove it from the system. Note that β is the global death rate, which we fix to be constant in time, whereas β/N (t) is the death rate for each particle. Since we only consider the steady state statistics in this paper, there is no distinction between fixing the global β or the individual particle's β/N (t) to be constant in time.
We shall use the following units for length [x] = d and time [t] = ζ/k. In these units, the particles' diameter is always 1 in simulation units. ζ/k is the typical timescale for two particles (separated by a small distance d) to separate. In our simulation units this timescale is again equal to 1. All the results below will be presented in simulation units.

III. RESULTS
We initialize all simulations from a single particle located at the centre of the box. We then let the system multiply and fill the entire box until a steady state is reached (see Supplementary Movies 1). Fig. 4(a) shows the total number of particles N (t) as a function of time t for one simulation run. Initially the total number of particles is one: N (0) = 1, and then N (t) grows exponentially until the steady state is reached at around t 50. At steady state, the total number of particles fluctuate around some average value. Some snapshots of the steady state configurations are shown in the insets of Fig. 5 steady state each particle has an average contact number roughly equal to or more than z max = 6, which corresponds to the maximal packing in two-dimension. At this packing, the particle does not divide anymore since the division rate α i is zero. Division is thus induced only if the particle dies (with finite death rate β/N (t)), at which point, the contact number drops below z max and the division rate becomes non-zero again.
We also measure the average velocity of the particles in the x-direction: The average velocityv x (t) as a function of time t for one simulation run with skewness value /L = 0.1 is shown in Fig. 4(b). It shows a clear bias of net velocity in the negative direction, indicating a macroscopic current induced by a combination of local division and apoptosis, and the asymmetric shape of the channel.
To make further analysis, we then take the time average and ensemble average ofv where t ss is the time it takes for the system to reach the steady state (in our simulations t ss ∼ 50). τ > 0 is the time length over which we take measurements (in our simulations, we take τ = 3t ss ). The angle bracket l e f t w a l l r ig h t w a ll 1 2 3 4 Figure 6. Net current of particles in the x-direction is mainly caused by division of the corner particles. By calculating the wall forces on the daughter particle, the average velocity vx ( ) as a function of can be computed. Green lines are parallel to each other, likewise red lines are parallel to each other.
. . . e in Eq. (6) indicates ensemble averaging over different simulation runs with the same set of parameter values. Fig. 5 shows the plot of v x as a function of skewness . Each data point is taken from an ensemble average of 100 independent simulation runs. In the symmetric case = 0.5L, the average x-velocity is zero as expected. For < 0.5L, we observe a net current in the negative x-direction, and for > 0.5L, we observe a net current in the positive x-direction.
It is worth noting that in the case of active extensile and contractile fluids such as bacterial suspensions and actomyosin network inside the cells [10], a macroscopic current can be observed even along a symmetric channel. In these examples, the particles are actually elongated and behave as 'force dipoles', which act on the surrounding fluid. A spontaneous symmetry breaking in the orientations of the particles can give rise to a macroscopic flow inside a symmetric channel in either direction (spontaneous flow transition [7][8][9]16]). In our case the particles are spherical and non-swimmers, thus the particles do not have any orientation (except when they divide into a dumbbell momentarily [17]). Thus we do not expect to observe any spontaneous flow transition in the symmetric channel i.e. = 0.5L. Now we will argue that the net current of particles in the x-direction, which we observe in Fig. 5, is mainly caused by particles' division at the corner walls. Let's consider the geometry of a corner particle in Fig. 6. Here, θ and φ are the slopes of the two adjacent walls. From Fig. 1, we can compute: Let's divide the region around the corner particle into four sectors: γ = 1, 2, 3, and 4 as shown in Fig. 6. Now suppose that this corner particle divides in the next time step. Then the probability that we place the centre of mass of the daughter particle inside sector γ is proportional to the angle of that sector γ. In particular, the probability that the daughter particle is placed inside γ = 1 or γ = 4 is (π − θ − φ)/2π and the probability that the daughter particle is placed inside γ = 2 or γ = 3 is (θ + φ)/2π. Now suppose that the daughter particle is placed inside sector γ = 4. Then the newly created particle will not experience any force from the walls. However if the daughter particle is placed inside sector γ = 2, the particle will then experience a wall force from the left wall. The x-component of this force is F w sin θ, where F w is the magnitude of the wall force. Similarly, if the daughter particle is placed inside sector γ = 3, the x-component of the force coming from the right wall is −F w sin φ. Finally, if the daughter particle is placed inside sector γ = 1, the particle will then experience forces from both the left and the right wall. The x-component of this force is: Therefore, on average, for every division of the corner particle, the walls will exert a force in the x-direction with magnitude: where we have substituted Eqns. (7)(8) in the last line above. Now we can assume that the average velocity of all the particles is proportional to the wall forces exerted on the system: v x ∝ F x,w (in other words, the tissue responds like a Newtonian fluid). We can then estimate the average velocity in the x-direction as a function of the skewness : Fig . 5 shows the comparison between the theoretical estimate (orange line) with the simulation data. In Fig. 5, the proportionality constant in Eq. (11) is obtained from the best fit value.
Although the simple theoretical model in Eq. (11) captures the direction of the current correctly, the curvature of v x ( ) does not fit perfectly to the simulation data, especially at extreme skewness values: 0.1L and 0.9L. This can be due to several factors. Firstly, the motions of the particles at high density are highly correlated. The force due to particles' divisions in the bulk can propagate to the walls [18]. Secondly, from Fig. 4(b), we also see occasional spikes of rather large negative net velocity. This is due to sudden jamming and unjamming of particles near the narrow gap between the top and the bottom wall, which were not accounted for in the theory Eq. (11).
We also measured the area fraction of the system in the steady state, which is defined to be the total area of the particles divided by the area available for the particles: where d is the particles' diameter and N is the average number of particles in the steady state. Fig. 7(a) shows the steady state area fraction ϕ as a function of division rate α 0 for fixed global death rate β = 2. As we can see from the figure, ϕ increases then saturates with increasing division rate as expected. Fig. 7(b) shows the steady state area fraction ϕ as a function of global death rate β for a fixed division rate α 0 = 20. Here, ϕ does not actually vary much with β, except for a very large value of β > ∼ 50 [dashed line in Fig. 7(b)]. In this regime β > ∼ 50, the death rate is so fast compared to α 0 so that the number of particles goes to zero.
From Fig. 7, we also note that the area fraction in the steady state appears to be larger than one: ϕ 1.8 for a large range of parameters. This is because a soft spring potential is used to model the interactions between the particles. This potential allows the particles to overlap, which mimic the mechanics of deformable particles such as foams [19]. Thus the diameter d in our model does not necessarily correspond to the diameter of the real cells (which are also deformable). The total overlaps between the particles also give rise to an internal pressure, which corresponds to the homeostatic pressure, observed in biological tissues [20,21].
Finally Fig. 8 shows the dependence of the average xvelocity v x on the division rate α 0 and the global death rate β. As we can see from the figure, v x is affected more strongly by the death rate β/ N than the division rate α 0 . As we argue in Fig. 6, the average velocity is mostly controlled by the division of the corner particles. At steady state, the corner particles have an average contact number close to z max = 6. Accordingly, the division rate of the particles is close to zero, see Eq. (4). Division is only triggered when one of the neighbouring particles dies with rate β/N (t), which then reduces the number of contacts. Therefore β sets the limiting rate for the  division of the corner particles.

IV. CONCLUSIONS
We have shown that local division and apoptosis events are sufficient to give rise to a macroscopic current (or work), without any self-propulsion mechanism. We have found that the maximum average velocity, obtained in our ratchet geometry is around 0.003 in simulation units. To convert this to physical units, we note that the typical cell size is around d 100 µm and the apoptosis rate is around β/ N 10 −6 s −1 [22]. We can then set the units of length and time in our simulations to be 100 µm and hour respectively. This translates to a maximum velocity of around 0.3 µm/hour, which is significantly smaller than that of bacterial ratchet. Thus although local division and apoptosis are sufficient to give macroscopic work, they are much less efficient than rectified self-propelled particles.
One way to improve the efficiency of our computational model is to consider the growth phase of the tissue rather than the steady state behaviour. In this paper we use a periodic boundary condition in the x-direction and the measurements are always taken in the steady state. However, at steady state, the division rate of the corner particles (which are responsible for the net current) is limited by the death rate β/ N . On the other hand, when we consider the growth phase of the tissue, where the total number of particles N (t) is still growing exponentially, the division rates of the corner particles are much faster and therefore we expect the net velocity to be much higher. To test this, we need to perform simulations in a sufficiently long channel without periodic boundary conditions and let the tissue grow into the empty space of this extremely long channel (see Fig. 9). We then expect the tissue to grow asymmetrically due to the ratchet effect, explained in this paper.
It might be interesting to replicate our ratchet geometry in real biological tissues (or growing bacterial colonies) to detect any net current such as the one depicted in Fig. 9. However, some cells such as epithelial cells are also motile [23] (in addition to local division and apoptosis events), which makes it difficult to disentangle the physics of local division and apoptosis from self-propulsion. Other type of cells such as most tumour cells [24] and biofilms [25], for example, are much less motile and may be more suitable to test our theory.
Finally, our results can have potential applications in tissue engineering, for example, in treating achilles ten-don injuries. One might engineer a ratchet-shaped scaffold starting from either end of the injured tendon, directed towards the middle, to speed up the healing process.

AUTHOR CONTRIBUTIONS
Both authors contributed equally to this work.

CONFLICTS OF INTEREST
There are no conflicts to declare.