Information transmission in a cell monolayer: A numerical study

Motivated by the spatiotemporal waves of the MAPK/ERK activity crucial for long-range communication in regenerating tissues, we investigated stochastic homoclinic fronts propagating in channels constituted by directly interacting cells. We assessed the efficiency of long-range communication through these channels in terms of the rate of information transmission. We identified stochastic phenomena that reduce this rate, most notably front propagation failures, new front spawnings, and variability in the front velocity. A tradeoff between the frequencies of propagation failures and new front spawning determines the optimal channel width (which geometrically prescribes the length of the wave front). We found that the optimal frequency of initiating new waves is determined by a tradeoff between the input information rate (higher for more frequently initiated wave fronts) and fidelity of information transmission (lower for more frequently initiated wave fronts).


Introduction
In living organisms, cells communicate through a variety of mechanisms, including chemical and mechanical signals.Long-range communication within a tissue may result from local communication between neighboring cells.This is the case of spatiotemporal MAPK/ERK activity waves, originating from the wound edge [1,2] or from leader cells [3], and involving a mechanochemical feedback loop that coordinates collective migration of epithelial cells [3,4].
Excitable ERK activity waves have been shown to control the rate of scale regeneration in zebrafish [5].These studies highlight the capability of waves to propagate recurrently across successive cell layers despite the inherent discreteness and heterogeneity of the communication medium.
From a dynamical systems perspective, a traveling front may be either an interface between regions in space that are in different equilibria (heteroclinic traveling waves in bistable systems) or an excitation that locally departs and then returns to a unique equilibrium (homoclinic traveling waves in monostable excitable systems) [6].The heteroclinic traveling waves are formed robustly at the interface between two different regions and as such are resilient to random perturbations [7].However, the passage of a heteroclinic wave irreversibly changes the state of the reactor that consequently cannot be re-used to support propagation of a subsequent front.In contrast, a homoclinic traveling wave is an out-of-equilibrium "stripe" flanked on both sides with the reactor in the equilibrium state.Although homoclinic traveling waves may be sent recurrently and in this way convey complex messages to spatially distant locations, they are fragile in the presence of stochastic fluctuations.Experimentally, the spatiotemporal ERK signal propagation has been observed to be distorted by random bursts of ERK activity [2,8], which occasionally give rise to spontaneous waves [9].
Here, we sought to quantitatively assess the capacity of a discrete, excitable medium to transmit information encoded within a train of activity waves.Specifically, we investigated stochastic homoclinic fronts propagating in narrow channels formed by directly interacting cells, and determined the efficiency of long-range communication through these channels in terms of the rate of information transmission, also referred to as bitrate.This metric quantifies the amount of information that can be transmitted through a communication channel in a unit of time.We identified several types of stochastic phenomena that reduce the fidelity and thus the rate of information transmission, among which the front propagation failure, new front spawning, and variability in the front velocity were the most impactful.We demonstrated that a tradeoff between the frequencies of front propagation failure and new front spawning determines the optimal channel width, enabling the fronts to reach the greatest distance and maximizing the rate of information transmission.We investigated the system's ability to relay periodic sequences of fronts as well as transmit binary-encoded information.Binary-encoded information is encoded by specifying a predetermined list of equally spaced time slots and deciding whether a front is sent or not sent at each of these time slots.We determined the time interval between the time slots that maximizes the information transmission rate.This optimal time interval or, equivalently, the optimal frequency of time slots results from a tradeoff between the input information rate (higher for on average more frequently initiated fronts) and fidelity of information transmission (lower for on average more frequently initiated fronts).Finally, our exploration of the model parameter space revealed that efficient long-distance information transmission is achievable only if the refractory time is several times longer than the neighbor-to-neighbor activation time.

Model definition
In our simulated monolayer, each cell assumes one of four distinct states: quiescent (Q), excited (E), inducing (I), or refractory (R), and inter-state transitions follow a predetermined . When activated by an I neighbor, a Q cell assumes the E state.Subsequent progression from the E to the I state enables the cell to induce activity in its Q neighbors.After the I cell assumes the R state, it loses the ability to activate its Q neighbors and becomes insensitive to activation by neighboring I cells.Finally, reverting to the Q state restores the cell's responsiveness to activation by an I neighbor.
Appropriate time scales of the residence in the E, I, and R states are crucial to enable wave propagation.In the context of the MAPK/ERK pathway, the first time scale is associated with signal reception and the signal transduction cascade (through the EGFR/SOS/RAS complex, RAF, and MEK), culminating in ERK phosphorylation; the second time scale is determined by the time required by phospho-ERK to trigger cell contraction that may result in activation of EGFR in neighboring cells [3]; and the third, longest, time scale is related to the refractory period of the signaling cascade, which, due to inhibitory multisite phosphorylation of SOS (by phospho-ERK), is at least partially insensitive to incoming signals.
In the model, the multi-step signal transduction within the MAPK/ERK pathway is reflected by the assumption that the E, I, and R states comprise multiple sub-states (nE, nI, nR, respectively).Consequently, the residence times in the E, I, and R states (E, I, and R, respectively) follow Erlang distributions, with their modes shifted away from zero and variances smaller than those of exponential distributions with identical means.the rates of transitions Q→ E1 (1/act), I1 → I2, and I2 → R (both nI/I) are all equal.Thus, the probability that an I cell will activate a given neighboring Q cell before transitioning to R is 0.5 + 0.5×0.5 = 0.75.This value is greater than both the site and the bond percolation thresholds on the triangular lattice, equal 0.5 and ~0.347, respectively [10], permitting front propagation [11].

Front speed
A propagating front consists of active cells (i.e., cells in either state E or state I) located in its head, followed by a thick block of R cells that prevent backward front propagation.We initiated fronts by setting the states of the cells in the first layer to I1 and observed the propagation of activity toward the other end of the reactor (Fig 1B).
In the model, the front is deemed to move one step forward once the next cell layer is activated and progresses through all E sub-states to become I.On average, a forward step, in which the front advances by just one cell layer, takes the time act/⟨nneigh_I⟩ + E, where ⟨nneigh_I⟩ is the average number of I cells in contact with a single Q cell.In a deterministic model, where the time spent by cells in each state is fixed, the front forms a straight line, and each Q cell at the front head has exactly 2 inducing neighbors.Thus, the propagation speed is vdeterministic = 1/(act/2 + E), which yields 1 layer/(4.5min) ≈ 0.22 layer/min for the nominal parameter values.
In the stochastic model, however, the front edge becomes rough and the cells to be activated have, on average, more than two I neighbors, which increases the propagation speed.We observed that for the nominal parameter values, the average front velocity ⟨v⟩ rises with the channel width W from vdeterministic for W = 1 to the asymptotic value vasymptotic ≈ 0.28 layer/min for W ≳ 10 (panel A in S1 Fig).

Transit time and its variance
The expected time in which the front travels the whole channel length L is ⟨transit⟩ ≈ L/⟨v⟩.For sufficiently long channels, the distribution of transit is nearly Gaussian with variance transit critically affects the fidelity of information transmission by determining the precision with which the moment of a front initiation can be inferred from the time at which it reached the end of the channel.

Propagation failure and front spawning
In a stochastic model, a traveling front is subjected to random events disrupting its propagation.In narrow and moderately wide channels, we observed two types of disruptive events: propagation failure and new front spawning (Fig 2A).
Propagation fails if all I cells progress to the R state before exciting any neighboring Q cells (S1 Video).New fronts are spawned when a cell remains in the I state long enough for one of its neighboring cells to recover from R to Q.Such a neighboring cell may get activated and become a source of a new front or fronts (S2 Video).The new front(s) can propagate backward or forward.When a backward-propagating front encounters a forwardpropagating front, they collide and usually annihilate (Fig 2B , S3 Video) or, rarely, give rise to another front.For broad reactors (>20 cells wide), fronts may propagate in directions not necessarily parallel to the channel longer axis (S2 Fig, S4 Video), which leads to a chaotic front pattern characteristic for the Greenberg-Hastings model (first defined by [12] and later recast as a stochastic model and studied, e.g., by [13]).This 2-dimensional effect is not observed in narrow channels, in which the front tail (the block of cells in the R state at the front's rear side) is longer than the channel width.
We found that the propensity (probability per one cell layer) of a propagation failure event decreases with the channel width W as where afail < 0 and Wfail are coefficients that depend on model parameters (Fig 2C).The exponential dependence on W results from the fact that the number of cells in either the E or the I state is proportional to the channel width, and for a front to disappear all of them have to simultaneously progress to R without exciting new cells.On the contrary, the propensity of new front spawning event increases linearly with channel width as with aspawn and Wspawn dependent on the model parameters.This is because spawning may be triggered by any cell across the channel width (Fig 2C).
We noticed that in narrow channels (W ≤ 6), a single backward front is usually spawned (Fig 2D

Optimal channel width
In the case of a narrow channel transmitting a series of fronts, the impact of the two types of disruptive events is similar.Propagation failure eliminates the front, whereas a backwardspawned front collides with and annihilates the subsequent (forward-propagating) front.In both scenarios, the total number of fronts received at the other side of the channel is reduced by 1.In narrow channels, a spawning event usually generates a single backward front.
Thus, the probability that a front in a series passes through the channel depends on the total propensity of disruptive events As fail decreases exponentially and spawn increases linearly with the channel width, there is an optimal width at which tot is minimized.This optimal width maximizes the probability of uninterrupted front propagation through the channel.For the nominal parameter values, the optimal channel width is Wopt = 6 (Fig 2C ), and for this width the average range of uninterrupted front propagation is about 8000 cell layers.As we will see later, the optimum is robust to moderate changes in the values of model parameters.We thus chose the value Wopt = 6 as the default channel width for most of the paper.Since uninterrupted front propagation is critical for information transmission, we expected that the information transmission rate would be highest when W is close to Wopt.

Refractory time
The tail of a front consists of a block of refractory cells, which hinder propagation of the next front.A cell requires, on average, Tcycle = act/⟨nneigh_I⟩ + E + I + R = 66.5 min to make a full cycle from Q through E, I, and R back to Q (assuming default parameter values and ⟨nneigh_I⟩ = 2).In the deterministic model, during this period, the channel is fully blocked; only fronts initiated (sent) at a rate smaller than 1/Tcycle are transmitted.In the stochastic model, Tcycle defines the time scale of the inter-front interval, below which the propensity of disruptive events is significantly elevated, as one can see in the kymographs in Fig 3A.
To investigate how often fronts can be sent and still be reliably transmitted, we performed simulations in which pairs of fronts were sent at various inter-front intervals (Fig 3B).As expected, for intervals shorter than Tcycle, the second front cannot be initialized successfully (initiation failure probability is 50% for the interval of 61.5 min), as the cells at the beginning of the channel have not recovered to Q yet.Additionally, even if a front is successfully initiated, it is much more likely to be eliminated later during the propagation.The proximity of the previous front also markedly increases the probability of front spawning.This is likely due to some R cells remaining after the passage of the first front, which recover and become Q only after the passage of the second front.These cells may get activated within the block of R cells in the second front's tail, seeding a new front at the second front's rear side.
The described effects were found to be most significant for inter-front intervals in the range of 50-130 min, i.e., around the value of Tcycle + cycle ≈ 96.6 min, where 2 is the variance of the time of the full cycle Q → E → I → R → Q.The time Tcycle + cycle =: TR can be considered the effective refractory time.For inter-front intervals longer than 130 min, the propensities of disruptive events approach the single-front values (shown in Fig 2C).For intervals shorter than 50 min, the propensities of disruptive events remain high, but the overall probability of their occurrence becomes low, because most fronts are not initiated successfully.

Propagation of periodically initiated fronts
The refractory time forces subsequent fronts to not travel too closely, and consequently the average interval between fronts reaching the end of the channel has a minimum (with respect to the frequency of initiating new fronts), which we denote Ttrans-min(L).For a short channel of length L = 30, the effective refractory time TR is a good approximation of Ttrans-min .Overall, we found that in the limit of short channel lengths, the empirically estimated refractory time TR sets an upper limit on the transmittable front frequency, while in longer channels, the maximum transmission frequency is lower, and lower initiation frequencies enable propagation of more coherent fronts trains.

Numerical results
To assess the bitrate, we used a simple binary protocol in which sequences of equiprobable binary digits {0, 1} are directly translated into sequences of fronts (1 → front initiated, 0 → no front initiated) and sent through the channel at regular time slots, with inter-slot interval Tslot.
We performed simulations and registered the times at which the fronts reached the end of the reactor.For each time slot tslot, we computed the expected arrival time texpected = tslot + ⟨transit⟩ = tslot + L/⟨v⟩ and selected the front that arrived closest to texpected.This could be a successfully transmitted front sent in the proper slot, a front sent in another slot (if in the proper slot a front was not initiated or during its propagation it was lost), or a front spawned in a disruptive event.We recorded the difference Δt = tarrival − texpected between the selected front's arrival time tarrival and texpected.Using all slots in the sequence, we estimated the mutual information (MIslot) between collected Δt values and the {0, 1} digits (see Methods for details).Finally, we calculated the information transmission rate r = MIslot/Tslot.The results for various channel lengths and inter-slot intervals are presented in Fig 4A .As one can observe, information transmission is highest for moderate values of Tslot.For long Tslot, the fraction of successfully transmitted information (equal to MIslot measured in bits, as each binary digit carries one bit of information) is determined by the distant-front dynamics and is thus roughly independent of Tslot (Fig 4B ).For large values of Tslot we may define two regimes: (1) The free-front regime, in which the interval between slots is at least twice longer than the transit time.In this regime, fronts spawned backward cannot collide with subsequent fronts, because they reach the beginning of the channel and disappear before the subsequent front is initiated.Consequently, information transmission is determined by the single-front propagation failure propensity.This regime is characteristic of short channels or very long inter-slot intervals.
(2) The distant-front regime, in which fronts are frequent enough to be annihilated by backward fronts spawned by their predecessors, yet still maintain sufficient distance to avoid direct interaction (for default parameters and L = 1000 this means Tslot > 200 min, see Fig 3C).In this regime, the propensity of disruptive events is still equal to that observed for single fronts, but information transmission is limited by the total disruptive event propensity, as both propagation failure and backward front spawning lead to the extinction of one front.
For the considered set of parameters, information transmission efficiency (MIslot) is substantially higher in the free-front than in the distant-front regime, because the annihilation of forward-propagating fronts by spawned backward-propagating fronts (Fig 2D ) is the main limiting factor in the distant-front regime.In both regimes, MIslot is nearly independent of Tslot, and thus the bitrate (equal to MIslot/Tslot) decreases as 1/Tslot.On the other hand, for short inter-slot intervals the bitrate is limited by strong interactions between fronts (increasing propensity of disruptive events) and transit time dispersion comparable to Tslot.Thus, for each channel length, there is an optimum inter-slot interval Topt for which the bitrate is the highest.
For the nominal parameter values and channel length L = 30, the optimum is located at Topt ≈ TR ≈ 96.6 min, and the maximum bitrate is ~0.5 bit/h.Unsurprisingly, the bitrate decreases with an increasing channel length L regardless of Tslot due to the accumulation of disruptive effects (Figs 4A, 4C).The optimal interval Topt increases with L, and in the investigated range of channel lengths, the increment is roughly proportional to √L (Fig 4D ), which can be attributed to the growth in the transit time variance.
The maximum information rate depends on the channel width (Fig 4E).As expected, the highest bitrate is observed for W = 6, for which the total disruptive event propensity is lowest (as shown Fig 2C).This maximum is lower but relatively more prominent for longer channels, when the impact of disruptive events on information transmission is strongest.
The dependence of the optimal inter-slot interval on the channel width is weak (Fig 4F ) and consistent with the changes in the maximum bitrate.

Semi-analytical predictions based on phenomenological analysis
There are two major phenomena that set bounds on the information transmission rate: disruptive events and transit time stochasticity.If the intervals between slots are long, the major limiting factor is the possibility of front extinction, due to either propagation failure or a collision with a spawned backward-propagating front.In this case, the amount of transmitted information per slot can be obtained from the confusion matrix with the formula MIslot = 1 − 1/2× [(p + 1) log2 (p + 1) − p log2 p], where p is the probability of front extinction (see Methods for details).To determine the value of p, we used the propensities fail and Thus far, the bitrate was estimated based on the arrival times of individual fronts.We show and discuss the inference based on two consecutive fronts, which yields higher bitrate estimates for short channels (L = 30), in S1 Text.

Disruptive events in single fronts
The maximum front frequency and information transmission rate depend on model We may notice that afail, and thus fail, do not depend on E and R, as those parameters have no influence on whether a cell will be activated.The failure propensity increases with act and decreases with I, because an increase of the ratio act/I implies a lower probability that an I cell activates a Q neighbor before proceeding to R.
The coefficient aspawn, and thus spawn, depend on all four kinetic parameters.It increases with act, which can be explained as follows: high act implies low cell activation propensity, which renders some cells activated at the rear side of the front and mediating front spawning.
Importantly, the value of coefficient aspawn increases with E (the increase of which implies a higher chance that the cell becomes inducible in the front tail) and decreases with R (the increase of which implies a broader zone of R cells behind the front).Another consequence of the discussed dependence of aspawn on E and R, and independence of afail of these two kinetic coefficients, is that an increase of R and/or decrease of E shift the minimum of tot(W) to higher W values.
In summary, we showed that an increase of the R/E ratio (i.e., the ratio of the residence time in the refractory state to the residence time in the excited state) nearly exponentially reduces the disruptive event propensity tot, and linearly increases the optimal width Wopt of the channel.

Periodic fronts and the bitrate
Understanding how the model kinetic parameters influence the propensities of disruptive events allows us to analyze and interpret their impact on Ttrans-min and maximum bitrate (Fig 6).The dependence of Ttrans-min on the act and I is dictated by the dependence of fail and spawn on these parameters; as expected, the increase of fail or spawn leads to the increase of Ttrans-min.The dependence of Ttrans-min on E and R is more intriguing: Ttrans-min has maxima with respect to E and decreases for large E; Ttrans-min decreases also for small R.This is puzzling, as we know from Fig 5A that large E and small R imply large spawn, which could block front propagation.However, as we can see in kymographs in Figs 6B, 6C, for small R and large E we observe a distinct pattern of front spawning.In this mode, numerous backward and forward fronts are created, so the front density (especially for small R) at the end of the channel is higher than at its beginning.In this regime it is impossible to coordinate collective cell motion in the direction of the front source.also suggest that there should exist some optimal R, at which the bitrate is highest.In fact, this maximum is attained close to the nominal value of R, which is 60 min (Fig 6D).Larger values of R do not allow for frequent fronts, while smaller ones are associated with a spawning front propagation pattern (blocking any information transmission).For the same reason, information transmission is blocked for large values of E.The information transmission rate increases monotonically with decreasing E.This is because the decrease of E reduces spawn and to some extent the refractory time, and does not influence fail.

Kymographs in Fig 6B
Unsurprisingly, in the MAPK/ERK pathway, E is short, despite ERK activation being a multistep process.As shown in S5 Fig, the maximum bitrate grows monotonically with the number of sub-states (providing that the total time of all sub-states remains constant).In agreement with the influence of act on fail and spawn, the bitrate also grows monotonically with act, and, similarly as with E, act appears to be short for the MAPK/ERK pathway.
Finally, the bitrate attains its maximum close to the nominal value of I.

Discussion
ERK activation triggers cell contraction, which leads to activation of EGFR in neighboring cells [3].In cell collectives, this mechanochemical coupling coordinates the propagation of the wave of ERK activity and cell movement in opposite directions [9].As the MAPK/ERK cascade is inhibited behind the front, cell contraction induces subsequent ERK activation in cells directly ahead of the front, rather than behind it.The mechanochemical coupling was theoretically studied in 1-D and 2-D models by Boocock et al. [4,14].Their model, constrained with data obtained from experiments on MDCK cells but devoid of the details of signal transduction through the MAPK/ERK cascade, allowed them to determine the optimal wavelength and period for maximizing migration speed towards the tissue boundary.In our study, we investigated the processes that interfere with the stochastic propagation of activity waves.
In single cells, the MAPK/ERK cascade has been found capable of transmitting information between membrane (opto)receptors and ERK at bitrate exceeding 6 bits/hour [15].This high bitrate results from the narrowly distributed delay between receptor activation and ERK phosphorylation, enabling well-synchronized activation of subsequent cell layers and robust propagation of ERK activity waves.Canonically, trigger (homoclinic) traveling waves employ a positive feedback to propagate over long distances [16][17][18].In our model, the positive feedback (at the tissue level) arises when the inducing cell excites a quiescent cell, and then the quiescent cell becomes inducing itself.Because after the passage of a homoclinic wave the system returns to its single steady state, such waves can be initiated recurrently at desired time points; the same property entails that in the presence of stochastic fluctuations waves can vanish but also may arise spontaneously.
In our study, we characterized two types of disruptive events: front propagation failure and new front spawning.Propagation failure eliminates the front, whereas a single backwardspawned front collides and annihilates with a subsequent front in the series.Thus, when fronts are initiated recurrently, in both cases, the total number of fronts is reduced by 1.
Consequently, the probability that a front in a series of fronts passes through the channel decreases with the total disruptive event propensity tot = fail + spawn.Importantly, because fail decreases (exponentially), while spawn increases (linearly) with the channel width W, there is some optimal channel width Wopt, for which the probability of uninterrupted front propagation through the channel is the highest.This result is surprising because the first intuition could be that the reliability of stochastic signal transduction monotonically grows with the channel width.
Each front on its rear side has a layer of refractory cells, which block the propagation of other fronts close behind them.The refractory time TR = Tcycle + cycle, where Tcycle is the cell cycle time (from Q through E, I, and R back to Q) and cycle is the cell cycle standard deviation.When studying periodically initiated fronts, we observed that the time TR sets the lower limit on the time interval between fronts that can be transmitted through a short channel.In longer channels, the minimum average interval between fronts Ttrans-min is larger than TR, and increases with the channel length.
To numerically estimate the rate at which information can be transmitted, we employed the binary encoding protocol.We found that for the optimal channel width, the bitrate is 0.5 bit/hour for L = 30 and 0.2 bit/hour for L = 1000 for nominal model parameter values accordant with the timescales of processes implicated in the MAPK/ERK signaling cascade.There are two major phenomena that limit the information transmission rate: disruptive events and transit time stochasticity.For long intervals between fronts, the possibility of front extinction (due to disruptive events) is the limiting factor.In this case, the amount of information transmitted per slot was found to be MIslot = 1 − 1/2× [(p + 1) log2 (p + 1) − p log2 p], where p is the probability of front extinction, which is an increasing function of front slot frequency.Since the bitrate is the product of the slot frequency and MIslot, it attains its maximum for some optimal inter-slot interval.When the interval between front slots is comparable to or shorter than the standard deviation of the transit time (that increases with channel length), fronts reaching the end of the channel may be attributed to a wrong slot.Once this effect is taken into account, our phenomenological predictions become satisfactory for both short and long inter-slot intervals.As expected, the bitrate reaches its maximum for the channel width which minimizes the propensity of disruptive events.
Finally, we performed a sensitivity analysis to show that the ratio of the refractory to the excited state residence times (R/E) nearly exponentially reduces the total disruptive event propensity.Surprisingly, we found that Ttrans-min decreases for large E and small R, attaining high values when the disruptive event propensity is high.Kymographs indicate that for small R or large E there is a distinct chaotic front spawning regime.In this regime, multiple backward and forward fronts are created, so that the front density (especially for small R) at the end of the channel is higher than at its beginning.Moreover, the existence of multiple backward fronts excludes coordination of collective cell motion in the direction of the signaling source.We found that (for other parameters fixed) there exists an optimal R (close to the nominal R value of 60 min) associated with the highest bitrate.Larger R preclude frequent fronts, while smaller R lead to the chaotic front spawning regime blocking any information transmission.Information transmission is also blocked for large values of E, and bitrate increases monotonically with decreasing E.This is because a decrease of E reduces spawn and does not influence fail.Unsurprisingly, in the MAPK/ERK pathway, E is short despite ERK activation being a multistep process.
Overall, our analysis indicates that there is an optimal channel width and optimal slot frequency, allowing the most reliable front propagation and consequently highest bitrate.
Additionally, sensitivity analysis suggests that the MAPK/ERK pathway is "optimized" for efficient information transmission by traveling ERK activity waves.

Numerical simulations
Monte Carlo simulations were carried out according to the Gillespie algorithm using the code adapted from [19].The simulator code is available at https://github.com/kochanczyk/qeirq The code that can be used to reproduce all reported results is available at https://github.com/pszyc/visavis-seir .

Front arrival time
Activity (i.e., total number of cells in states E or I, NEI) in the last cell layer was counted as a function of time; local maxima along the temporal axis were found and filtered based on peak height to remove peaks closer than 20 or 30 min apart.

Disruptive event detection
Activity (i.e., the total number of cells in states E or I, NEI) was counted as a function of the distance x from the beginning of the channel and time to form a kymograph (x-t plot; Figs 3A, 6B, 6C).To determine front position(s), the kymograph was smoothed (with a Gaussian kernel along the x-axis and an exponential kernel along the t-axis); local maxima of NEI along the x-axis were found (thresholds on minimal peak height and minimal distance between peaks were used to discard short-lasting fluctuations).Fronts were tracked using LapTrack [20].Fates of fronts that did not reach the end of the channel were classified as 'annihilation', 'initiation failure', or 'propagation failure' depending on whether, respectively, the front had disappeared in the vicinity of another front propagating in the opposite direction, within a certain minimal distance from the channel beginning (1-5 cell layers depending on the channel width), or while propagating freely through the channel.To detect front spawning, track splits were recorded.Tracks shorter than 50 min that did not split into other tracks were discarded together with the split in which they were created.The remaining track splits were treated as front spawning and grouped into events: two splits were assigned to a single spawning event if they occurred not further than 20 cell layers apart in space and 150 min away in time ( Figs 2D, 2E).

Front speed and transit time variance
To determine the front speed and transit time variance (S1 Fig) , for each channel width, N = 30,000 simulations with a single front were performed in a reactor of length L = 300 (or L = 30 in the cases when, due to a very short average front range, too few fronts reached the distance of 300 to collect sufficient statistics).The average and the variance of the first front arrival time were used to compute the front speed and the transit time variance.

Information transmission ratesemi-analytical predictions
To predict the information transmission rate for distant fronts (Fig 4B ), we used the coefficients aspawn, afail, Wspawn, and Wfail obtained based on single-front simulations (Fig 2C ; Eqs (1-2)) to calculate spawn and fail.Then, we computed the probability that a front is eliminated  = exp (− × ( fail +  ×  spawn ×  backward )), where nbackward = 1.285 is the expected number of backward fronts spawned in a single event and information transmission rate as To include the effect of interaction between fronts (Figs 4G-4I, dotted line), we took into account the probabilities of disruptive events pfail = ppropagation failure + pinitiation failure and pspawn obtained from the numerical simulations with two fronts at different intervals (Fig 3B ), as a function of the initial interval.To account for the fact that in the binary protocol the initial interval between fronts can be any multiplicity of the inter-slot interval Tslot (interval of length k × Tslot has probability 2 −k ), we averaged the probabilities pfail and pspawn using the formula where pevent is either pfail or pspawn.Then, we approximated the probability that a front is eliminated with the formula where nbackward = 1.285 was taken from the single-front simulations and  was calculated as Based on p, we calculated MIslot and bitrate using Eqs (6-7).Note that, unlike before, in this approach, we used probabilities of disruptive events rather than propensities, which required running the two-front simulations for each channel length separately (analogous to Fig 3B, in which results only data for L = 300 are presented).This was necessary, as the probabilities do not scale linearly with the channel length, due to the fact that the disruptive events most likely take place near the beginning of the channel.
To include additionally the effect of the transit time variance (Figs 4G-4I, dashed line), we assumed that the transit time has a normal distribution with  2 =  2 transit =  2 0 L and computed the probability of a front being attributed to a wrong slot (i.e., arriving further than Tslot/2 from the expected arrival time) pinaccurate = 2 − 2 (Tslot / (20√L)), where  is the c.d.f. of the N(0, 1) distribution.We modified the elimination probability p' = p + (1 − p) pinaccurate and computed the probability that a front was mistakenly detected in a slot in which no front was sent pfake = (1 − p) pinaccurate − 1/4 (1 − p) 2 pinaccurate 2 (the latter term prevents double counting of cases in which fronts from both the previous and the next slot are attributed to the considered slot).We created the confusion matrix: where c00 = (1 − pfake)/2, c01 = pfake/2, c10 = p'(1 − pfake)/2, c11 = [(1 − p') + p'pfake]/2, and computed MIslot using the standard formula: Figures  to the probability that a corresponding number of backward and forward fronts was spawned in a single event; the total area of all disks is proportional to the front spawning propensity, which is different for each channel width.
E Propensity of spawning one or more fronts in a single localized disruptive event as a function of channel width.
Data for panels C-E was gathered from 30,000 simulations for each channel width, channel length was fixed at L = 300.In all panels, channel width W = 6.Data for panel B from 3000 simulations for each data point; data for panels C-E from 30 simulations of 500 fronts for each interval.

C, D
Maximum bitrate (C) and optimal inter-slot interval (D) as a function of (the square root of) the channel length.Channel width W = 6.

E, F
Maximal bitrate (E) and optimal inter-slot interval (F) as a function of the channel width for various channel lengths.

G-I
Prediction of bitrate taking into account probability of front extinction due to propagation failure or collision with a backward front (dotted line), and the chance that a front was attributed to a wrong slot due to transit time stochasticity (dashed lines)see Methods for details.Colored lines as in panel A.
Data for panels A-B and G-I from N = 100 simulations with 500 front slots for each data point.Data for panels C-F was computed by searching the maximum of a graph computed as in panel A.

S1 Text
Alternative decoding.In the main text, the bitrate was estimated based on the arrival times of individual fronts.One may expect, however, that an inference based on several consecutive fronts may yield higher bitrate estimates.The use of auxiliary fronts could for example facilitate discerning whether a given received front originated from a true sent front or rather was born out of a disruptive event.As a step towards exploration of higher-order decoding schemes, we recomputed the bitrate estimate based on the arrival time of the front closest to texpected and the arrival time of the preceding front.For longer reactors, the bitrate increased only slightly (see figure below), but for shorter ones, the estimates are markedly higher, especially for shorter inter-slot intervals.For L = 30, wherein ⟨transit⟩ ≈ Topt, a secondary local bitrate maximum appears close to Topt/2.In the considered case, decoding based on two arrival times helps distinguish whether a front was not received when expected due to it not being sent or because the channel was blocked by a front sent in the previous slot.This demonstrates that it is possible to increase the bitrate despite a high likelihood of losing the subsequent front if two fronts are sent in consecutive slots.
In a simulated fully-confluent monolayer, cells are immobile agents arranged on a 2-D triangular lattice of length L and width W (Fig 1B).For nominal parameter values (Fig 1C), ).The value of  2 0 decreases with the channel width (panel B in S1 Fig).The variance  2 For a long channel, L = 1000, we observe that Ttrans-min is about 200 min and is achieved when fronts are sent every 160 min, while a more frequent front initiation results in a slightly longer average front arrival interval.In long channels, a considerable percentage of fronts sent at periods < TR is eliminated shortly after initiation (Fig3D).For fronts sent every 150 or 180 min, the propensity of front elimination is initially low, but increases slightly with the distance from the initiation site (Fig3D).This is because for such intervals, disruptive events are relatively rare, but due to fluctuating velocities some fronts draw closer together, and thus the likelihood of disruptive events increases.Consequently, in a channel of length L = 300, fronts sent every 150 or 180 min arrive with time spans distributed around the initial period (Fig 3E, top row), whereas fronts sent with periods < TR arrive with time spans distributed more broadly around Ttrans-min regardless of the initial pulse frequency (Fig 3E, bottom row) Fig 4B.By replacing the distant-front event probabilities with estimates for finite interval between fronts from Fig 3B we obtained satisfactory predictions of the bitrate for Tslot > 120 min for L = 100 and L = 300 and for Tslot > 200 min for L = 1000 (Figs 4G-4I, dotted line).
parameters.In Fig 5A we analyze how the kinetic model parameters influence the propensities of disruptive events, fail and spawn, given in Eqs (1-2).For the nominal parameter values, Wfail ≈ 0 and Wspawn ≈ 1, and numerical analysis indicates that these coefficients remain in the range (−1, 2) for the considered range of model parameters (S4 Fig), and thus for the nominal channel width, W = 6, has a modest effect on fail and spawn.Therefore, crucial for understanding the system's behavior are the changes of coefficients afail and aspawn.It is important to notice that afail changes approximately linearly with the kinetic parameters, while aspawn scales exponentially; as a consequence, both fail and spawn exhibit exponential dependence on E, I, R, and act.
Coefficient aspawn decreases with nE, nI, and nR, because the increase of the number of intermediate states renders the distribution of times of transitions E → I, I → R, and R → Q narrower, resulting in less stochastic front propagation.Consequently, the minimum value of tot is an increasing function of E, decreasing function of R, and decreasing function of nE, nI, and nR (Fig 5B).

(
obtained from data shown in Fig 2D) and  is the probability that a backward front collides with a next front before reaching the channel beginning.We used  = 1 to obtain the distantfront value shown in Fig 4B.Based on the probability of front extinction p we constructed the confusion matrix with entries TP = (1 − p)/2, TN = 1/2, FP = 0, FN = p/2 and computed MIslot as

Figure 2 .ABD
Figure 2. Disruptive events associated with front propagation.

Figure 3 .BCDE
Figure 3. Occurrence of disruptive events as a function of the interval between fronts.

Figure 4 .
Figure 4. Rate of information transmission.

S5Fig.
Bitrate dependence on the number of E, I, and R substates.The encircled dots in each panel correspond to the nominal parameter set given in Fig 1C.

Bitrate estimation based on
more than one front arrival time.Bitrate estimation is based either on the timing of the front that arrived closest to the expected arrival time or of the closest front and the preceding one.Figure analogous to Fig 4A.