In brief, the model is multi-compartment, composed of a two-dimensional grid that represents the inner cellular membrane connected to a three-dimensional reservoir, that represents the cytoplasm. Several molecule types are coexisting and interacting in the framework (G-actin, actin crosslinkers (ACs), and myosin), and can transit between the grid and the reservoir. The total number of molecules for each type (cortex + cytoplasm) remains constant throughout the simulation. For simplicity, diffusion in the cortex and other processes such as rotations are not allowed.

The code is written in Julia language [37], taking advantage of the following libraries: StatsPlots, Dates, KernelDensity, LightGraphs, FFTW, Plots, Random, Printf, Statistics, GLM. Code is available upon request.

Numerical simulations are performed using a Coarse Grained framework developed in house. The model is solved numerically using a Monte-Carlo computational algorithm. The dynamics for each event is accepted or rejected based on the Metropolis method [38]. The probability distribution of the different micro-states of the systems is evaluated statistically using a semi-grand canonical ensemble approach.

In brief, the model works as follows: 1) G-actin molecules in the cortex are allowed to polymerize directionally to form filaments (F-actin). 2) F-actin is then stabilized and linked by ACs, forming meshworks of parallel, anti-parallel and perpendicular filaments. 3) Myosin can bind to F-actin, and produce mechanical force. 4) When bound to two anti-parallel F-actin, the mechanical action of myosin results in tension or stress. 5) When the stress supported by a given F-actin is higher than its attachment to the membrane, the F-actin is released from the cortex to the cytoplasm, and part of its tension is redistributed to the rest of the network. Details and parameters for each of these processes are explained below (a more exhaustive and detailed explanation of the model can be found in ref. [31]).

At every time step, a location and orientation (horizontal or vertical) in the grid is chosen. If the location is empty, an actin monomer (G-actin) can transition from the cytoplasm to that location in the cortex based on the value P + = e − μ 1 μ ₁ represents a chemical potential with the following form: μ 1 = μ 1,0 − γ 1 ∗ l o g N 1 N 1,0 (1) where μ _1,0 is a reference potential, N _1,0 is the total number of actin molecules in the system, N ₁ is the number of actin molecules in the grid at a given time point. Parameter γ ₁ is a modifier of the slope of the function.

The potential function μ ₁ sets the probability of transition of a given molecule from the reservoir to the cortex. The change is accepted or rejected based on the Metropolis method. In the most rigorous definition of the chemical potential in thermodynamics, the units are energy units per molecule), but in the context of semi-grand canonical ensembles, the chemical potential is defined with no units.

To allow crossing of F-actin, each position in the grid can accept one G-actin molecule in vertical orientation and one in horizontal orientation. If insertion is successful, a link between G-actin and the cortex is established, with energy E ₁ (mimicking the effect of integrin molecules [19]). In addition, if another G-actin exists directly behind and with the same orientation, a link between the two molecules is established (with energy E ₀), and the filaments grows directionally.

On the contrary, if the location selected is already occupied, removal of the G-actin is evaluated using again the Metropolis method ad the value of P ⁻ = e ^−E , being E the sum of the energy of the links that affect the molecule. Directional de-polymerization is implemented by spontaneous release of the bond of the last G-actin of a F-actin (therefore, P − = e − E 1 in this case).

ACs can occupy positions between pixels of the grid, so they can join filaments. Similarly to G-actin, we start by selecting a random location in the grid. If location is between two parallel, anti-parallel or perpendicular filaments, insertion is evaluated based on P + = e − μ 2 and the Metropolis method, with μ ₂ being a chemical potential for ACs: μ 2 = μ 2,0 − γ 2 ∗ l o g N 2 N 2,0 (2)

If insertion is accepted, links of energy E ₂ are established between the AC and the two F-actin.

If the location already has an AC, removal is evaluated based on the value P − = e − 2 E 2 , and the Metropolis method. To study the properties of the actomyosin network, two classes of ACs are used: ACs that only link parallel F-actin, and ACs that can also link perpendicular F-Actin. A scheme of how the two mechanisms of action is shown in Figure 1 (right panel).

Myosin molecules work in our model similarly to ACs, with their own chemical potential μ ₃ and their link energy E ₃. μ 3 = μ 3,0 − γ 3 ∗ l o g N 3 N 3,0 (3)

When linked to anti-parallel F-actin, myosin adds a certain amount of stress every time step. When the stress accumulated by a filament is higher than the energy attaching it to the cortex, it will be removed, and a fraction of its stress is communicated to the filaments that were connected to it. This way, the stress propagates through the network and it is progressively disassembled. This mechanism is at the core of the oscillatory behavior that we will analyze in the Results section.

Changes in shape of a cell will affect the shape of its neighbors. To implement numerically this spatial coupling, we run simultaneously two numerical simulations corresponding to two neighboring cells (cell i and cell j), and make the value of a given parameter dependent on the total stress tension (TS, the sum of all tension in all filaments of the cell). For the coupling of the potential functions, we define the following function: μ 1 , j = μ 1 , j u n c o u p l e d 1 + α β γ T S i γ + β γ (4)

Another approach is to modulate the stress in 1 cell by the stress of the neighbour cell (if 1 cells compressed, the other experiences a relaxation in the network tension): to do that, we define the following function: T S j = T S j u n c o u p l e d 1 + α β γ T S i γ + β γ (5) in both previous equations α affects the maximum intensity of the coupling, β affects the inflection point and γ affects the slope. This way, the value of the parameter uncoupled in cell j is recalculated based on the value in the neighbor cell i, using the function that goes from 1 (no effect of neighbour) to α (full effect of neighbour). In all simulations, we use values of α = 5, γ = 5. These values have been chosen to provide an intermediate level of coupling (smaller values show very little effect, higher values suppress the oscillations in both cells). β is set as the half mean value of the total tension, which for the parameters used moves from 0 (minimum of oscillation), to around 10 (maximum total tension at the peak of the oscillation), therefore β = 5 also.

Spatial and temporal calibration of this type of coarse grained models using experimental data is challenging because each interaction in the model often combines many biochemical steps (by design, the model is a very simplified version of the real system, and many molecular interactions are condensed as one simple step). The chemical potential is used as a tool that we use to maintaining a dynamic equilibrium between the molecules in the grid and in the reservoir, therefore it does not have a direct comparison with any experimental value (parameters of the chemical potential are chosen to maintain intermediate levels of cortex occupancy). Finally, bond energies are defined in terms of multiples of RT (being R the universal gas constant). Bonds mediated by ATP hydrolysis are on the order of 30 kJ/mol (RT = 2.47 kJ/mol). The power of Myosin in breaking the network is calibrated based on the time it takes in experiments to disassemble the actomyosin cortex (based on our own reference [30]).

To understand the main features and properties of the oscillations in both main components of the actomyosin network (actin and myosin), we start with a characterization of the dependence of the dynamics on the model parameters, to define which aspects of the molecular interaction between actin, myosin and ACs affect the amplitude, period or robustness of the oscillations.

Figure 2A focuses on the effect of the reference chemical potential μ _1,0, which regulates the affinity of the G-actin to transit from the cortex to the grid (therefore, it effectively affects the concentration of G-actin in the grid versus the cytoplasm). We can see that oscillations appear for an intermediate value of μ _1,0 (low values result in very small amounts of G-actin unable to polymerize, while high values result in too high affinity of G-actin towards the cortex). In this situation, new G-actin enters the grid before the network is fully disassembled, and the system is driven out of the periodic regime towards a dynamic steady state situation (filaments are continuously being assembled and disassembled, but no periodicity arises as a global property of the system.

This is also illustrated in the left panel, where we plot the power of the highest peak of the corresponding Fourier transform of the amount of actin in the cortex overtime, computed for more values of μ _1,0 (right panel). The maximum periodicity occurs for values of μ _1,0 around −3.

Figure 2B focuses on the effect of the link energy between G-actin and the inner cell membrane E ₁ in the period and amplitude of the oscillations. We can see that oscillations are more clear for an intermediate value of E ₁ (low values of E ₁ result in very unstable filaments; high values of E ₁ result in filaments too attached, so myosin cannot produce the sufficient tension to release the filament, so disassembly occurs as a more stochastic and less periodic event). The representation of the maximum power of the Fourier spectrum, shows that the optimal values for periodicity occur around E ₁ = 7 RT, but it is more or less maintained for very high and very low values of the same parameter.

Finally, Figure 2C represents the effect of the power provided by myosin (W ₃), which regulates the energy that a given myosin motor produces per unit of time, once attached to two anti-parallel F-actin. Changes in these parameters seem to strongly affect the period of oscillations, with much longer and higher amplitude bursts as we move towards less power of myosin. This makes sense, since it takes longer to reach the threshold tension to detach a given filament. The longer it takes, more time for the filaments to grow in size and to fill the grid (therefore, amplitude is also increased). Again, optimal periodicity also occurs at intermediate values (see right panel for a maximum peak of the Fourier transform).

A more exhaustive analysis of the robustness of the oscillatory dynamics for each parameter of the model is presented in Table 1, where we include a brief statement explaining the biological meaning for each parameter, and the minimum and maximum values for clear and sustained oscillations. Model parameters used to produce this figure and the rest of the figures in the manuscript are liste in Table 2.

Since oscillations of the actomyosin cortex are in essence the assembly and disassemble of a interconnected network, we can use this analogy to study the properties of such network from a graph theory perspective. In this approach, we can define F-actin as the nodes of the network, and the links are the molecules capable of attaching to F-actin, i.e., ACs and myosin. This way, a poorly connected network represents a situation where F-actin molecules are linked to a small number of other filaments, while a highly connected network represents a situation where F-actin is linked to multiple other polymers.

To quantify this connectivity, we use the shortest path length L (also called geodesic distance, i.e., the minimum number of links required to connects two nodes in a network), calculated using the Dijkstra algorithm [40]. The average shortest path length < L > corresponds to this value averaged for all shortest paths from each node to each node in the network. For simplicity, we use the maximum connected component size transformed in a simple graph, without weights. To minimize the effect of networks of different sizes, we divide by the average shortest path length < L 0 > for an equivalent network of the same number of nodes and links, but where connections are established randomly. Therefore, high values of < L / L 0 > represent a highly ordered network, while values of < L / L 0 > around 1 represent a network structure similar to random.

Since our inner cell membrane is simplified as a simple square grid of discrete pixels, this will strongly affect the architecture of the network. To study different configurations and different situations of the graph, we first analyze a condition of networks composed of only parallel and anti-parallel F-actin (see Methods and Figure 1). This situation illustrates the structure and properties of bundles linked only by Fascin and similar ACs [41].

Figure 3A represents this first scenario. The left panel is a plot comparing the value of < L / L 0 > as a function of time (purple line), with the oscillations in the amount of actin in the cortex (green line). We can see that the value of < L / L 0 > increases and decreases in time, concomitant with the cortex occupancy, suggesting that as the network grows, it differs more from a random-like connectivity (since only parallel and anti-parallel filaments can be connected, two filaments that are placed far away from each other will have a very long L).

If we compute the degree distribution at the point of minimum (central panel) and maximum (right panel) of the oscillation, we see that in both situations, the maximum probability is slightly above two (since F-actin can only bind in this scenario to parallel or anti-parallel neighbors, most of the filaments are only linked to two other F-actin, and only a few filaments are linked to a higher number of F-Actin). The inserts in each panel correspond to spatial representations of the graph at these two time points. The network develops with few connections per node towards a highly regular configuration at its maximum size (occupancy ≈ 80 , % < L / L 0 > ≈ 4 ) (the ring form occurs because the boundaries in our model are defined as periodic).

Figure 3B represents the scenario of F-actin being able to link also perpendicularly (see Methods and Figure 1). This situation represents the properties of actomyosin networks where ACs are also able to link perpendicular F-actin polymers, mimicking the role of ACs such as the αActinin family [42]. The value of < L / L 0 > (purple line, left panel) again oscillates in time, but now the average shortest distance decreases as the network size increases i.e., as the network grows, it resembles more a random network. This behavior corresponds to the opposite scenario of Figure 3A. The explanation for this difference is the following: when the network is beginning to form, the most probable connections are between parallel and anti-parallel filaments (perpendicular filaments can only be joined by their intersection point). These initial connections result in an increase in < L / L 0 > , similar to the case in Figure 3A. However, as soon as connection between perpendicular filaments appear, < L / L 0 > falls drastically to a value near one, meaning that the average shortest path length is much more similar to what would be expected in a random graph.

This different connectivity impacts how tension and energy is being redistributed in the network, which in the end affects the amplitude of the oscillations (green line, left panel). The addition of perpendicular links results in larger oscillations that become more asymmetric in their shape, compared to smaller and more sinusoidal oscillations in parallel networks (Figure 3A, left panel). This asymmetry in the shape of the oscillation occurs because the higher connectivity makes the tension travel faster and further across the network, resulting in a more sudden and organized disassembly.

The central panel in Figure 3B shows the degree distribution at the minimum occupancy point, with an average degree slightly above 2 and small clusters being formed (similar to panel A). The right panel in Figure 3B corresponds to the maximum point of occupancy, with a much broader distribution (average degree around 6) and a much larger network structure (more connections, see insert for a graphical representation of the network in graph form), resulting in a connectivity very close to random.

The actomyosin cortex is an active adaptive material, capable of a wide range of dynamic behaviours, from bulk pulse contractions [26] to anisotropic polarized flows [43–45]. These traveling waves of actomyosin play pivotal roles in cell division, cell motility, and morphogenesis. Wound healing experiments in Drosophila show perhaps the clearest example of these cortical flows, with a single pulse of actomyosin filaments assembling from cell to cell toward the margin of the wound [46]. Another key example of actomyosin waves occurs during contraction and protrusion in Dictyostelium discoideum [47].

In our previous publication, we have demonstrated how our minimal model is able to reproduce the bulk oscillations of the actomyosin responsible for the periodic contraction of cells. In this section, we take advantage of the fact that our model has been constructed as spatiotemporal to study the conditions that will result in the arising of polarized actomyosin flows.

Figure 4A represents three snapshots of the system for conditions of bulk oscillations. We can see how the cortex transits from full (t ₁) to less crowded (t ₂) to almost full again (t ₃). A time lapse movie of the system for these conditions is presented as Supplementary Movie S1. On the other hand, if we reduce the value of E ₁ (Figure 4B), instead of bulk oscillations, the system starts to exhibit traveling waves that swipe the system directionally. In this regime, the total occupancy levels do not change so drastically along the simulation, but still we observe oscillations in local levels if we focus on a small region of the system. A time lapse movie of the system for these conditions is presented as Supplementary Movie S2.

To further illustrate the difference between these two dynamic regimes, we generated space-time plots (Figure 4C), where we pile the intensity levels for the same pixels at different time points of a given simulation. The left panel represents a situation of bulk oscillations, where all the actomyosin is disassembled from the cortex simultaneously. From this plot, it can be seen that the periodicity is highly regular since very early in the simulation, and the collapse of the network is almost instantaneous (illustrated by a red horizontal line).

The right panel corresponds to the situation of low E ₁. Now, instead of full sudden collapse of the network, the system transits into a regime where disassemble occurs initially disordered in small clusters that, as time progresses, become more organized in the form of spatio-temporal waves.

Next, we compare these two potential scenarios with the dynamics observed in experiments performed in the basal region of Drosophila follicle cells, as well as in amnioserosa cells during dorsal closure (see Methods). In brief, to analyze myosin dynamics in the basal side of Drosophila follicle cells, we performed live imaging of egg chambers expressing the in vivo marker for myosin, Sqh:mCherry. To analyze myosin dynamics in the apical side of amniserosa cells, we use Drosophila embryos stably expressing zipperCPTI002907 to label myosin.

To do this, we recorded time lapse movies and monitor the changes in intensity levels and the spatial redistribution of the actomyosin cortex. Figure 4D illustrate the dynamics in the experiments in the form of a space time plot (the full movie for myosin oscillations in follicle cells is presented as Supplementary Movie S3, the full movie for myosin oscillations in amnioserosa cells is presented as Supplementary Movie S4). To monitor the oscillatory dynamics in the experiments, we compute the levels of myosin instead of F-actin, because myosin staining allows us to distinguish better single cells (as seen in our previous publications [30], the dynamics of both actin and myosin are very similar).

Left panel in Figure 4D corresponds to a single cell quantification of the follicle tissue (Supplementary Movie S3). The changes in intensity levels show a dynamic that initially resembles wave-like oscillations (marked with oblique lines as a guide to the eye) but later in time, the dynamics is more bulk-like (there is no clear direction in the assembly and disassembly). On the other hand, in amnioserosa, (right panel, the waves can be clearly distinguished during the 20 min of the experiment (see Supplementary Movie S4). Interestingly, although the speed of the waves in both tissues is around the same order (slope of the lines), the waves are very different in shape (myosin waves in amnioserosa are much thinner, i.e., disassembly occurs almost immediately after assembly, while in basal follicle cell oscillations, the disassemble takes longer, suggesting that tension accumulates more slowly). Based on the model analysis, we hypothesize that both experimental tissues represent different in regions of the phase space: basal follicle cell oscillations are more similar to oscillations of high attachment (Figure 4A), while amnioserosa oscillations resemble a situation of low energy attachment. In addition, we have to take into account the inherent experimental differences between both systems (see Methods), and that neighboring cells are also oscillating and probably affecting the dynamics of their neighbours. This effect of the coupling between cells in close proximity is studied in the next section.

There is vast evidence that confirms that actomyosin oscillations are cell autonomous, rather that arising as a dynamic property of the system. On the other hand, since cells are not isolated, but coupled mechanically as part of a larger tissue, it is safe to assume that the dynamics of a given cell will be influenced by the dynamics of their neighbours. Since the integrity of epithelial tissues such as the amnioserosa or follicle cells has to be maintained during morphogenesis, the expansion (or contraction) of a cell in a particular direction has to be coupled with a complementary contraction (or expansion) of its neighbours. In this direction, a clear experimental evidence of the effect of neighbors in the dynamics of the oscillation has been described during wound healing, where the process of cell expansion is simply a result of the absence of tension by the actomyosin, combined with tension being generated in neighboring cells [48].

Previous studies focused on the inter-relation between the changes in cell shape at the tissue level show that the axis of maximum contraction tend to organize in their orientation, showing for the first time a supra-cellular organization during dorsal closure [8]. This is consistent with actomyosin waves that can expand several cell diameters, but authors did not focus on this characterization.

A more detailed studied of the spatiotemporal coordination of the pulses during the invagination of the Drosophila ventral furrow, establishing a link between the myosin intensity and the change in area of the cell [49]. Here, authors find that in wild type conditions, the tissue is not governed by force balance, meaning that the reduction in area of a given cell is not affected by peaks in myosin levels in neighboring cells.

To investigate if the mechanical interaction between neighboring cells in the tissue affects in some way the dynamics of the actomyosin cortex, we use again experiments in Drosophila follicle cells (see Methods). We perform cell tracking of a set of cells in the tissue and measure the intensity levels of myosin overtime. A time lapse movie of the process is presented as Supplementary Movie S5. Next, for each cell, the oscillations are fitted to a sinusoidal function to extract parameters of the dynamics, such as amplitude and phase of the oscillation. Finally, these two parameters are compared between neighbors to identify potential coupling or the existence of supra-cellular organization.

Results of the quantification are summarized in Figures 5A–C. Panel A represents a snapshot of the time lapse movie, with the labeled cells part of the analysis. Figure 5B corresponds to the oscillation frequency of all the cells analyzed, plotted in levels from blue (low) to red (high frequency). Our data suggests that cells coordinate in the form of arrays with similar frequency of oscillations (illustrated by orange lines). This supra-cellular organization in the form of long lines of cells is similar to the arrays presented in ref [8], and suggest that, a similar organization to the compression and expansion is also taking place at the level of the frequency of oscillation. Interestingly, the frequency values (lower panel) show a distribution that resembles a bimodal shape, but a more extensive analysis (more cells, more movies) will be required to conclude this with mode confidence.

Supra-cellular actomyosin waves should result in spatial coordination in the phase of the oscillations between neighboring cells. This quantification is presented in Figure 5C, where the phase relative to cell 0 is represented in a color-bar (restricted to cells with a similar frequency, since comparison of phase is only informative when performed between cells with similar dynamic characteristics). Interestingly, cells closer to our reference cell have a similar phase, while cells further in the tissue oscillate in anti-phase with the reference cell. These results points towards the existence of close coupling in terms of actomyosin levels, suggesting that, despite the lack of coordination in area changes (see ref [49]), the actomyosin levels seem to be coupled between neighbours.

Next, we go back to our model to study potential scenarios for inter-cellular actomyosin coupling that can reproduce the experimental observations. To do that, we generated two parallel numerical simulations of 2 cells simultaneously, and tested different scenarios for the coupling. Details of how the coupling is implemented are described in the Methods section. The first scenario corresponds to coupling in via chemical potential. The rationale is that, in a tissue, the contraction of cell i should result in expansion of cell j. Therefore, an increase in area may correspond to a reduction in local concentration (represented in our model by the reference chemical potential μ _0,i ). The actin levels of the 2 cells in this situation are represented in Figure 5E, where we mark with black and red rectangles the regions where the 2 cells oscillate in phase and anti-phase, correspondingly. Comparing this situation of μ _0,i coupling with Figure 5D (no coupling), we cannot see real differences in the phase (rather than an decrease in the frequency of both cells). Increasing the dependence between both cells seems to suppress the oscillations, rather that increasing the synchronization. In addition, to quantify how this coupling we affect the synchronization between the 2 cells, we compute the correlation matrix for each pair of simulations and compare the correlation terms. For the case where no coupling is present, the correlation obtained is -0.02 ± 0.08 (blue bars in right column of Figure 5D). When we add the coupling term to the calculation of the chemical potential, the value is similarly low (-0.09 ± 0.17), implying that there is no coupling or that the periods of phase and anti-phase are very similar and cancel each other.

Next, we performed simulations where the coupling is set to occur at the level of the tension. The rationale behind this scenario is the following: since a compression of a cell results in expansion of the neighbor, this expansion may correspond to a release in the local tension of the network (due to local breakup of links between network and inner cell membrane, for instance). The resulting oscillations in this scenario are represented in Figure 5F. Compared to the control simulations with no coupling (Figure 5D), there seems to appear more regions where locking between the two phases occur. Interestingly, the synchronization is more marked in phase (black rectangles) that in anti-phase (red rectangles), which agrees with the experimental observation that cells in close proximity seem to be oscillating in phase (Figure 5C, blue region). Correlation analysis when the systems are coupled through their total stress results in a value of 0.42. ± 0.12, meaning that this approach is causing a more important effect than the previous two (blue bars in right column of Figure 5D). Again, higher intensity in the coupling seems to suppress the oscillation, rather than increase the synchronization.