Mechanical oscillations in biological tissues as a result of delayed rest-length changes

Mechanical oscillations in biological tissues as a result of delayed rest-length changes Jos´e J. Mu˜noz1, M´onica Dingle1, Manuel Wenzel2 October 30, 2018 [email protected], www.lacan.upc.edu/munoz 1. Dept. Mathematics, Universitat Polit`ecnica de Catalunya-Barcelona Tech 2. Karslruhe Institute of Technology Abstract Oscillatory behaviour in tissue biology is ubiquitous and may be observed in the form of either chemical or mechanical signals. We here present and solve a mechanical cell model which exhibits oscillations emerging from delayed viscoelastic rheological laws. These include a time delay between the mechanical response and the rest-length changes, which evolve proportionally to the delayed cell deformation and use a remodelling rate parameter. We show that different regimes (no oscillatory response, sustained oscillations, and unstable oscillations) are obtained for different values of the delay or the remodelling rate. The results are analytically demonstrated in a one-dimensional problem with one and two cells which are represented by simple line elements. Oscillations of the cell deformations are obtained whenever different delays co-exist, or the delay is size-dependent. We also extend our results to a multicellular two-dimensional vertex model that includes the same rheological law, and which inherits the presence of critical values of the delay or remodelling rate. We numerically show that indeed the size-dependent rest-length changes induce oscillations in the cell shapes and areas. 1 Introduction Oscillatory responses in tissues are observed in numerous phenomena: gene regulatory signalling [1], cardiac pulsation [2], shear waves in cardiac tissue [3], chemical synapses in neural loops [4] or pattern formation during embryo segmentation [5, 6], to name a few. The suggested physical causes behind the observed oscillations are also varied: cardiac electrical impulses [2], cycling gene expression [6], inertial terms [3], or delays produced by distance in signalling and responses [5, 4]. In turn, different methods have been used to model these systems: activator-repressor coupling in lattices [1], second order differential equations [3], control differential equations with feedback [2], diffusion-reaction 1 equations [7, 8], or delayed contractile forces [9, 10]. Despite these studies, it is yet unclear how these ubiquitous oscillatory responses are originated in in vitro [11] or in vivo epithelia [9]. In this study we concentrate our attention on displacement oscillations observed in embryonic tissues [10], and model them through a delayed response. Previous use of a delay for modelling oscillators can be found for instance in [5, 9], and has been justified by the presence of a distance between the sender and the receiver of the biochemical signal, or due to the time needed to process such signal. We will not investigate the causes behind this delay, but rather analyse its consequences on the oscillatory behaviour and stability of the cell deformation. We show that the delay strongly affects the stability, which in turn may help to explain embryonic development. In embryo morphogenesis, area oscillations appear precisely prior to drastic shape changes of the organism, and have been suggested as a potential mechanism originating ventral invagination [12], dorsal closure [9, 13] or symmetry breaking [14] in Drosophila fly. The dynamics of these cell area fluctuations have a strong resemblance to wave-like oscillations of a fluid surface (see movies in the supplemental material of [9] and Figure 1 in [10], which illustrate the oscillations at the apical surfaces). However, in cellular systems, inertial terms are negligible, which renders the governing mechanical partial differential equations as elliptical or parabolic, with no apparent oscillatory solution. It is thus necessary to search for other sources that may trigger such fluctuations. Recently, these oscillations have been correlated with contraction pulses and myosin concentration pulsations [15, 16, 17]. Moreover, the time lag between mechanical displacements and concentration oscillations has been experimentally measured, and is approximately equal to one quarter of the wave period (∼230 s) [10]. Motivated by these observations and experimental measurements, we here resort to a minimal system that includes elastic mechanical equilibrium and a delay in the cell active response. The latter is represented by an evolution equation of the cell rest-length that gives rise to a delay differential equation (DDE). DDEs have been extensively used in control systems [18] and have been recently explored in cellular oscillators [5]. In the context of embryogenesis, contractile forces with delays have been also employed in [9]. However, the stability and oscillatory conditions have not been analysed for such cellular systems. The coupling between cortex turnover, its contractility, mechanical equilibrium and the measured oscillations have been studied for a two cell system in [19], while the associated myosin transport and cell area changes are modelled in [10]. In contrast to these works, we decide to omit the solution of transport diffusion-reaction equations. In our model, the delay in the evolution equation is attributed to the retarded bio-chemical signalling, and simulates the coupling of the myosin concentration and the cortex contractility. By directly inserting the delay between the cell mechanical state and the active response, we manage to retrieve analytical conditions for the stability as a function of this delay and the cell material properties. 2 2 Active rheological model with delay 2.1 Single element analysis We here model the cell with a line element that has a simple linear elastic law. The element is schematised in Figure 1a, and has a stress given by σ=kεe, with εe=`−La displacement based strain measure, and `and Lthe lengths at the deformed and stress free configuration. In order to mimic a viscoelastic response, we resort to the following evolution law for the rest-length L[20]: ˙ L=γ(εe(t)−εc).(1) Contractility parameter εccorresponds to a material homeostatic value at which Lis constant [21]. It has been shown that this law, in conjunction with the linear elastic law mimics a Maxwell model with viscosity η, where the characteristic time k/η is here given by the remodelling rate γ[20]. Similar models based on rest-length changes have been employed for embryonic tissues in order to represent their ability to adapt and remodel [22, 23]. For simplicity in our subsequent analysis, we assume that εc= 0, although the validity of the methodology and conclusions remain unchanged without this assumption. We further hypothesise that the rest-length adapts with a given delay τ, so that equation (1) adopts the following expression, ˙ L=γ(l(t−τ)−L(t−τ)).(2) (a) (b) Figure 1: (a) Schematic of model for 1 cell. (b) Model with 2 cells. We now impose a constant stretch u=l−l0, with an initial rest-length L(0) = L0=l0, τ ≤t≤0. The resulting delay differential equation (DDE) and initial conditions reads: Find L(t),∀t > 0,that satisfies ˙ L=γ(l−L(t−τ)) with constants l > l0, γ > 0, τ > 0 and initial conditions L(t) = l0,−τ≤t≤0. (3) The stability of the DDE can be analysed resorting to the solution of the characteristic equation [18, 24, 25]. It is known that there is a critical value τoscil of the delay for which the rest-length L(t) oscillates, and another value 3 τstab above which the solution is unstable, giving oscillations with increasing amplitude. We will show that in our case these limiting values depend on the remodelling rate γ. We will recall the main results and methodology, which will be recalled in the analysis of the 2 element system shown in Figure 1b. . We will assume that the solution can be written as a linear combination of exponential forms L(t)m= (L0−l0)emt +l0. By inserting this form in the DDE in (3), with m∈C, we obtain the characteristic equation m+γe−mτ = 0.(4) By setting m=α+iβ, we arrive to the following system of equations: α+γe−ατ cos(βτ) = 0 (5) β−γe−ατ sin(βτ)=0.(6) Equation (5) reveals that the system admits infinite solutions (α, β)k, k ∈Z for each pair of values γ, τ > 0. Since the DDE is linear, the solutions are given by the general expression: L(t) = l0+ (L0−l0) +∞ X k=−∞ emkt. Stable solutions are obtained for α < 0, and oscillatory solutions whenever β6= 0. It can be demonstrated that the following condition holds for oscillatory solutions [18], τγ > 1 e,(7) while unstable solutions are present whenever τγ > π 2. These results allow us to define oscillatory and stability limits given by, τoscil =1 eγ , τstab =π 2γ(8) Also, for α= 0, βτ =π 2±nπ, and thus, at the limit of stability, the solution is oscillatory with frequency π 2τ. Furthermore, for β= 0, we obtain α=−γe−ατ <0, and therefore non-oscillatory solutions are stable. We have solved the DDE in (3) numerically for some values of γand τ. The resulting rest-lengths are shown in Figure 2 and confirm the analytical results. Figure 3 shows the stability diagram, which indicates the regions with oscillatory and stable solutions in the (τ, γ)-plane. 4 0 20 40 t 0.8 1 1.2 L(t) =1.1* stab =3* oscil =oscil (a) 0.8 1 1.2 L -0.2 0 0.2 dL/dt =1.1* stab =3* oscil =oscil (b) Figure 2: One element system. (a) Time evolution of rest-length L(t) , with one example of non-oscillatory solution (thick dashed red line), a stable oscillatory solution (thick green continuous line), and a unstable oscillatory solution (thin blue continuous line). The three solutions correspond respectively to region I, II and III in Figure 3. (b) Phase diagram of L(t) for the three solutions . III II I 12345 1 2 3 4 5 Figure 3: Stability diagram for 1 element. Solutions in region I are nonoscillatory (τ < τoscil), in region II are oscillatory and stable (τoscil < τ < τstab), while in region III are oscillatory and unstable (τstab < τ). Thick lines correspond to functions τγ = 1/e (boundary regions I and II) and τγ =π/2 (boundary regions II and III). 5 2.2 Two element analysis We here analyse the two elements system in Figure 1b. Both ends of the system are fixed and only the node between the elements is allowed to move horizontally with a displacement u(t). The equilibrium condition between the nodes reads σ1=σ2, which for k1=k2=kis equivalent to the kinematic condition, l1−L1=l2−L2.(9) with liand Lithe the apparent and rest-length of element i= 1,2, respectively (see Figure 1b). The time evolution of the rest-lengths of elements 1 and 2 are given by, ˙ L1=γ1(l1(t−τ)−L1(t−τ)) ˙ L2=γ2(l2(t−τ)−L2(t−τ)) After using the fact that l1+l2=l0=const and the stress equilibrium in (9), the previous equations may be expressed as, ˙ L1=γ1 2(l0−L1(t−τ1)−L2(t−τ1)) (10) ˙ L2=γ1 2(l0−L1(t−τ2)−L2(t−τ2)) (11) This is a set of coupled DDEs. We will analyse its stability and oscillatory response for the simpler case γ1=γ2=γ. As before, by suggesting solutions of the type Li(t) = li(0) + (Li(0) −li(0)emit,i= 1,2, we get similar but coupled stability conditions, m1+γ 2e−m1τ1+e(m2−m1)te−m2τ1= 0, m2+γ 2e−m1τ2e(m1−m2)t+e−m2τ2= 0. When searching for time-independent values of mi, we must have m1= m2=m, and thus we arrive to one single equation: m+γ 2e−mτ1+e−mτ2= 0. The oscillatory and stability conditions can be now studied from the system of non-linear equations α+γ 2e−ατ1cos(βτ1) + e−ατ2cos(βτ2)= 0,(12) β−γ 2e−ατ1sin(βτ1) + e−ατ2sin(βτ2)= 0.(13) For non-oscillatory solutions (β= 0), equation (12) implies that α < 0, so that non-oscillatory solutions are stable. When β→0, equations (12)-(13) turn 6 into α+γ 2e−ατ1+e−ατ2= 0,(14) 1−γ 2e−ατ1τ1+e−ατ2τ2= 0.(15) By setting p1=ατ1and p2=ατ2, and multiplying (15) by α, the oscillatory condition can be found by searching pairs (p1, p2) that satisfy e−p1+e−p2≥ − e−p1p1+e−p2p2(16) and setting τi=2pi γe−p1+e−p2, i = 1,2.(17) The stability condition can be determined by setting α= 0, in which case equations (12)-(13) read cos(βτ1) + cos(βτ2) = 0 (18) β=γ 2(sin(βτ1) + sin(βτ2)) (19) From the first equation we conclude that oscillatory solutions (β6= 0) imply that βτ1=βτ2+π, which inserted in (19) gives the stability condition: π τ1+τ2 ≥γsin πτ1 τ1+τ2 .(20) The stability diagram in Figure 4 indicates the regions with oscillatory and stable solutions for L1and L2for different values of γand in the (τ1, τ2)-plane. Due to the kinematic and mechanical coupling, the solution is equally oscillatory or stable for each one of the two elements, and has a larger instability region as γincreases. We also point out that even if the delay of one element satisfies τi> τstab, the solution may be stable if the other delay is sufficiently below τstab. Similarly, even if τi> τoscil, solution may be non-oscillatory when the delay of the adjacent element is below τoscil. It can be said, that the delay of one element may help to stabilise or damp the global solution. The direction along the axis τ1=τ2in each stability chart of Figure 4 corresponds to an horizontal line in the stability charts of Figure 3 for the values γ= 0.5 and γ= 1. Indeed, for τ1=τ2, equations (16) and (20) give the same oscillatory and stability condition for one single element in (8). This can be also deduced from the fact that when the delays in both elements are equal, we have that ˙ L1=˙ L2, and the sum of rest-lengths L=L1+L2evolves 7 III II I 1=2=0.5 0 5 1 0 2 4 2 (a) III II I 1=2=1 0 5 1 0 2 4 2 (b) Figure 4: Stability diagram of 2 element system for different values of γin the (τ1, τ2)-plane. Solutions of L1and L2in region I are non-oscillatory (τ < τoscil), in region II are oscillatory and stable (τoscil < τ < τstab), while in region III are oscillatory and unstable (τstab < τ). Contour plots have been obtained with numerical solution of DDE, while thick lines correspond to analytical solution in equations (16)-(17) and (20). as ˙ L=γ(l0−L(t−τ)), so that the same critical values of the delay apply for L1+L2. For τ1=τ2we also have that ˙u=d dt(l1−l1(0)) = 1 2 d dt(l0−(L1−L2)) = 0 and thus no oscillations appear in the displacement of the cell boundary. This fact, and the presence of cell area oscillations in experimental measures motivates the use of a size-dependent delay with the form τ(t) = λl(t).(21) We could not retrieve in this case the analytical critical values of λfor stable and oscillatory condition (DDE becomes non-linear). However, we have numerically tested the relation in equation (21), and confirmed that displacement oscillations are obtained for a size-dependent delay, with the same value of λfor the two cells, but with unequal initial sizes. The stability diagram in Figure 5 compares the region boundaries of one element (thick lines) and the one using a size-dependent delay. In the latter case, the values of τin the horizontal axis correspond to the value of λ. The initial lengths of each element are equal to l1(0) = l2(0) = 1.0. For these lengths, the stable region is slightly reduced in the (λ, γ)-plane when compared to the stable region in the (τ, γ)-plane. Figure 6 shows the time evolution of the rest-lengths and displacements for the case with constant but different delay in the two elements (left), and for a size-dependent delay (right). Both give a oscillatory displacement u(t), but for the values of τand λtested, the amplitude of the average oscillation diminishes. 8 III II I 12345 , 1 2 3 4 5 Figure 5: Stability diagram for 2 elements. Thick lines correspond to regions τ1=τ2(same lines as in Figure 3). Contour regions correspond to sizedependent delays τi=λli, obtained with li(0) = 1. 0 10 20 30 40 t -1 0 1 2 L1 L2 u (a) 0 10 20 30 40 t 0 1 2 3 L1 L2 u (b) Figure 6: Time evolution of length Land displacement ufor: (a) different delays τ1= 0.5, τ2= 5 and (b) size-dependent delays τi=λli, with λ= 2.5. The two solutions correspond to the stable oscillatory regions (II) in Figure 5. We note that the use of a size-dependent delay is motivated not only by the triggering of oscillations at the cell boundaries, but also by potential physical hypotheses. For instance, 1) biochemical signalling between the nucleus and the cortex may take longer on larger cells, 2) the effects of (de)polymerisation of actin may take longer to build up, or 3) a globally contractile cortex requires more time to form on larger domains. None of these explanations has been confirmed, but are possible explanations for justifying the more realistic response of a size-dependent delay. 9