Spatial heterogeneity in collective electrotaxis: continuum modelling and applications to optimal control

Collective electrotaxis is a phenomenon that occurs when a cellular collective, for example an epithelial monolayer, is subjected to an electric field. Biologically, it is well known that the velocity of migration during the collective electrotaxis of large epithelia exhibits significant spatial heterogeneity. In this work, we demonstrate that the heterogeneity of velocities in the electrotaxing epithelium can be accounted for by a continuum model of cue competition in different tissue regions. Having established a working model of competing migratory cues in the migrating epithelium, we develop and validate a reaction-convection-diffusion model that describes the movement of an epithelial monolayer as it undergoes electrotaxis. We use the model to predict how tissue size and geometry affect the collective migration of MDCK monolayers, and to propose several ways in which electric fields can be designed such that they give rise to a desired spatial pattern of collective migration. We conclude with two examples that demonstrate practical applications of the method in designing bespoke stimulation protocols.

Since the experimental data contain information on the bulk velocity after stimulation with the electric field is turned off, it follows that the observed decaying velocity profile can be used to immediately infer the decay rate, γ.From these considerations, we seek to fit an exponential of the form to the bulk velocity data in the 2.5 hours after stimulation with the electric field is turned off.
Here, t end is the time at which the electric field is switched off, and C is given by We perform a least-squares fit for the decay rate, γ, to give an estimate of γ = 1.765h −1 .Figure A shows that the fitted exponential decay curve is in excellent agreement with the experimental data after pulse stimulation.This motivates us to take the value of γ as fixed going forward.The benefit of directly estimating γ from data is that it reduces the dimensionality of the parameter space, thus simplifying inference of the remaining parameters.

S2 Understanding recoil in the top edge post-stimulation
Here, we turn attention to understanding the phenomenon that the velocity in the top edge of the tissue becomes negative for several hours after the electric field is switched off.This finding was originally reported by Wolf et al. [1] who refer to it as recoil.In fact, recoil was also mentioned by Wolf et al. for the velocity at the trailing edge, which exhibits a sharp drop below its equilibrium value before restoring.On this phenomenon, Wolf et al. proposed that the observed recoil might be due to a mechanical effect induced by the visco-elastic properties of the epithelium.Although such a mechanical effect may indeed be at play at the trailing edge, we reason that it is not sufficient to explain the sustained negative velocity at the top tissue edge, as the characteristic time scale for the elastic response of MDCK epithelia is on the order of 15 to 30 minutes [3].Therefore, the recoil as observed in the experiments of Wolf et al. [1] cannot be explained by visco-elastic mechanical effects alone in the top edge.We then ask the question: if mechanics are not responsible for creating recoil, might it be that cessation of the electric field signal itself could be inducing cells to move in the opposite direction of the electric field post-stimulation?

S2.1 Model and constitutive equations
So far, we have assumed that the effective signal inside the constituent cells making up the monolayer is always positive given a positive field, as we have considered s + eff as in Equation (17) in the main text.The effect of this is that the effective signal inside the cells can only give a migratory command in the direction of the electric field.To understand whether cessation of the electric field itself can provide a negative signal to the cells, we wonder if the model assumption that the effective signal is always in the direction of the field can be relaxed, i.e. by replacing s + eff by s eff in the active force for the top edge.Since we only consider the dynamics of the top edge in this case, we can considerably simplify the model at the top edge, since there is no need to compare the response of the top edge relative to that of the bulk.Without loss of generality, we can set δ = 1, and study a reduced system of equations describing the dynamics of the top edge of the tissue, )

S2.2 Bayesian parameter identification
As before, we keep the pre-fit value for γ constant, and seek to estimate the parameter values for α, τ e , and τ a using the same Bayesian inference procedure as in Section 3 of the main text.We find that the posterior means of the model parameters are given by 0.103h, 6.11h, and 57.37h variations in τ e result in much smaller overall changes in the solution.In particular, reducing the value of τ e leads to very small changes in the solution, with the rise of the velocity to its peak being nearly indistinguishable between different values of τ e in the limit as τ e → 0. We argue that, therefore, non-identifiability of the parameter τ e in this case can be easily understood given the fact that these variations in model solutions are much smaller than the variance present in the experimental data.On a practical level, these issues arise because data acquisition in the experiment is only in 10 minute intervals;xt the data are not temporally resolved enough to tell different values of τ e apart when there is only one time series.The inference results in Section 3 of the main text demonstrate that this issue does not occur in the presence of more data, but in the case of a single time series the temporal resolution becomes problematic.
In summary, we have shown that relaxing the assumption that the effective signal in the monolayer is always positive can explain the finding that the velocity in the top edge of the monolayer turns negative for several hours after the electric field is turned off.By carefully investigating model predictions with varying parameter values, we found that the undershoot is well-described by the maximum likelihood value of the parameters, which was difficult to identify in practice owing to the sloppiness in the model.Interestingly, such an undershoot is not observed any of the other tissue regions studied in this work, which begs the question for

Figure A :
Figure A: Experimental data of bulk velocity decay post-stimulation with an electric field pulse of 3 V/cm (black dots) together with least-squares fit of an exponential decay curve in the form C exp(−γt) (solid blue line).The least-squares solution gives an excellent fit to the data and can be used to estimate the decay rate, γ.

Figure C :
Figure B: Bayesian inference for the top edge velocity model in Equation (S3).Left panel: posterior 95% confidence interval for the top edge velocity model in green with experimental data (black scattered points).Right panel: posterior distributions for model parameters α, τ a , and τ e , respectively, from top to bottom.