Optimal Control of Collective Electrotaxis in Epithelial Monolayers
Simon F. Martina-Perez, Isaac B. Breinyn, Daniel J. Cohen, Ruth E. Baker
TL;DR
This paper introduces a framework to optimally control the movement of epithelial cell layers using electric fields, enabling precise steering of collective cell migration.
Contribution
The paper presents the first unified framework for optimal control of collective electrotaxis in epithelial monolayers.
Findings
An ordinary differential equation model was developed to predict monolayer velocity in response to electric fields.
Optimal electric field designs were derived to maximize distance, velocity, or maintain constant velocity.
The framework can be extended to other external cues for steering collective cell migration.
Abstract
Epithelial monolayers are some of the best-studied models for collective cell migration due to their abundance in multicellular systems and their tractability. Experimentally, the collective migration of epithelial monolayers can be robustly steered e.g. using electric fields, via a process termed electrotaxis. Theoretically, however, the question of how to design an electric field to achieve a desired spatiotemporal movement pattern is underexplored. In this work, we construct and calibrate an ordinary differential equation model to predict the average velocity of the centre of mass of a cellular monolayer in response to stimulation with an electric field. We use this model, in conjunction with optimal control theory, to derive physically realistic optimal electric field designs to achieve a variety of aims, including maximising the total distance travelled by the monolayer, maximising…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8- —http://dx.doi.org/10.13039/100000893Simons Foundation
- —http://dx.doi.org/10.13039/100000009Foundation for the National Institutes of Health
- —http://dx.doi.org/10.13039/100000001National Science Foundation
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsItalian Literature and Culture · Linguistic Studies and Language Acquisition
Introduction
Electrotaxis is the process by which eukaryotic cells establish a cell polarity and move directionally in the presence of an electric field (Cohen et al. 2014; Allen et al. 2013; Zajdel et al. 2020). Electrotaxis is strongest in cellular collectives due to E-cadherin mediated adhesions (Wolf et al. 2022; Li et al. 2013), and displays striking similarities between different cell types and electric field stimulation patterns. Previous studies have focused on using electrotaxis as a means to achieve particular patterns in collective movement (Cohen et al. 2014; Zajdel et al. 2020). For instance, electric fields can be programmed to give bespoke migration patterns (Zajdel et al. 2020), which are along the field lines in converging and diverging electric fields (Cohen et al. 2014). A common theme in all of these experiments is that the average direction of collective migration can be accurately controlled during electric field stimulation. The ease with which electric fields can establish directed migration in cellular collectives has led electrotaxis to be accepted as a robust method of steering collective migration (Leal et al. 2023). The theoretical interest in modelling collective electrotaxis then lies in the ability to answer the following question: Given a desired spatio-temporal movement pattern, how should the electric field be designed to achieve this?
The adoption of electrotaxis as a predictable driver of collective migration has been stalled by a poor understanding of the temporal dynamics of the collective speed during electric field stimulation. For example, Wolf et al. (2022) found that MDCK collectives display a strong, but transient response to the electric field: while the average speed in the direction of the field increases quickly after the electric field is turned on, the collective starts to slow down after approximately one hour, even though the electric field strength remains constant. In general, the current lack of understanding of how the speed of the cellular collective will evolve during stimulation with an electric field poses difficulties in designing bespoke stimulation strategies that can be used to achieve a predictable migration outcome. This difficulty places the question of predicting collective speed during electrotaxis among the most pressing issues in collective electrotaxis. The problem we seek to address in this work is therefore how the collective speed of MDCK monolayers during stimulation with a given electric field can be predicted, and how this knowledge could be exploited to design a stimulation strategy to produce a desired pattern of collective migration?
Both theoretically and experimentally, the question of how to control the temporal dynamics of the collective velocity has scarcely received attention in the literature, even though the fact that the temporal dynamics of electrotaxis are non-trivial is well-documented. Controlling the temporal dynamics of collective velocity remains challenging, since the subcellular processes governing the evolution of cell polarity and active force production in response to stimulation with an electric field remain poorly understood (Allen et al. 2013; Cohen et al. 2014; Wolf et al. 2022; Zajdel et al. 2020). At the subcellular level, electrophoresis of charged membrane components creates a cell polarity that activates intracellular signalling pathways that are largely shared with those found in chemotaxis (Allen et al. 2013; Song et al. 2021; Gao et al. 2015). Subsequent modelling efforts Ifunanya Nwogbaga et al. (2022) have been able to explain the directional movement of cells during electrotaxis, but the existing literature is not well positioned to explain the temporal dynamics of cell responses.
In this work, we are interested in constructing a deterministic, continuous-time model to predict the average velocity of the centre of mass of a cellular monolayer in response to stimulation with an electric field. We will consider only uniaxial electric fields, i.e., the field can be described solely by its field strength in the positive or negative direction. Under these assumptions, the goal of this work will be to derive a system of ODEs which predict the collective velocity and the active forces in the monolayer. A significant benefit of using ODE models is that there exists a completely developed theory on how to optimally control ODE solutions according to a given optimisation problem Lenhart and Workman (2007). In the context of collective electrotaxis, this means that a well-calibrated and reliable ODE model can be used to solve an optimisation problem, whereby the optimality condition is given by the desired temporal pattern of collective movement and the control is the electric field input.
The current scholarship on collective migration of epithelia provides a convenient starting point for the construction of such an ODE model – see Alert and Trepat (2020) for a review of continuum models of collective epithelial migration. For example, since the internal signaling pathways for collective electrotaxis are largely shared with chemotaxis, one can begin to model the internal cellular response to an electric field with a well-established continuum model for the response to a chemotactic signal. Erban and Othmer proposed an adaptation-excitation model for the chemotactic response of single cells, whereby cells are excited when initially exposed to an external signal, and slowly adapt over time Erban and Othmer (2005). In their work, signal transduction is interpreted as “having two input pathways, an excitatory one that stimulates the [response], and one [...] which, in turn, shuts off the response”. Similar models are frequently proposed to describe chemotaxis. The advantage of such a model is that it is given by a straightforward system of scalar ODEs that describe the effective signal transduced within a cell upon contact with a chemoattractant (Erban and Othmer 2005). In this work, we will assume that similar internal signaling dynamics can be used to model the cellular response to an applied electric field. By letting the strength of the active forces in the monolayer be described using the model of Erban and Othmer (2005), we can derive a closed system of ODEs to describe the monolayer speed. Using experimental data in conjunction with Bayesian inference, we can interrogate the validity and applicability of using such a system of ODEs to model collective electrotaxis in response to stimulation with an electric field. Put together, the model can be used to develop a framework for designing electric fields that lead to desired outcomes in the form of an optimal control experiment.
This work will be structured as follows. In Sect. 2 we present the experimental data of collective electrotaxis that will be used for calibration and validation of the ODE models. Then, in Sect. 3, we derive a simple ODE model that describes the evolution of the tissue velocity. We then use this system of ODEs to predict the optimal control policy when a number of targets are considered, and compare our results to the current standards for electric field stimulation used in the literature. Together, we offer a mathematical modelling framework to describe, understand, and optimise key cellular processes involved in the regulation of collective electrotaxis.
Methods and Data
In this section, we present the publicly available experimental data from Wolf et al. (2022). In brief, the experiments of Wolf et al. (2022) involved electrotaxis of 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {mm}^2$$\end{document} MDCK-II epithelial monolayers. In these experiments, a total of nine MDCK epithelial monolayers were grown to confluence and stimulated with an electric field of 3V/cm—see stimulation trace in Fig. 1. In each experimental replicate, 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times 5\hbox {mm}^2$$\end{document} square tissues were seeded onto 10 cm tissue-culture plastic dishes (CELLTREAT) using stencils of 250 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} m-thick silicone elastomer (Bisco HT6240) cut with a hobbyist razor-writer (Silhouette Cricut). Suspensions of MDCK cells were seeded into the stencil patterns at a volume of 10 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} L and density of 2.5 (±0.15) x 10^6^ cells/mL. Cells were given one hour to adhere before the dish was flooded with media and incubated for 18 h prior to stencil pull and insertion into the electric bioreactor. Tissue culture stencils were pulled approximately 1 h before imaging began. Immediately post-stencil pull, dishes were refilled with fresh culture media followed by direct integration with the bioreactor. The assembled bioreactor is then placed into the microscope, situated within a custom-built cage incubator maintained at 37^∘^C and perfused with fresh media continuously bubbled with 5% \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_2$$\end{document} using a peristaltic pump (Instech Laboratories) at a rate of 2.5 mL/h. Electrotaxis was induced by custom electro-bioreactors developed by the Cohen group at Princeton University, which are based on the SCHEEPDOG platform (Zajdel et al. 2020). These bioreactors provide a continuous flow of media and can be programmed to generate electrical fields with custom electric field densities.Fig. 1. Experimental data of MDCK monolayer electrotaxis with its electric field stimulation protocol. Left panel: phase field image of an MDCK epithelial monolayer indicating direction of the uniaxial electric field and bulk region where the mean velocity is measured. Right panel, top: electric field density during the experiment. Right panel, bottom: bulk velocity corresponding to the electric field trace in the top row
Velocity data from each of the experiments were acquired using particle image velocimetry (PIV) on phase contrast images, where the cell nuclei are used as the fiduciary markers. Local tissue velocities were computed every 10 min. The PIV data provide a velocity field at a spatial resolution of 32 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 32 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {m}^2$$\end{document} , which corresponds roughly to 3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 3 cell lengths. For a sample velocity heat map, see Fig. 1. For further details on the PIV data processing and implementation details in FIJI Schindelin et al. (2012) and MATLAB The MathWorks Inc. (2022), we refer the reader to Wolf et al. (2022). Experimental PIV data provides detailed spatiotemporal information on the velocities at a fine spatial resolution. Wolf et al. (2022) previously used this detailed spatial data to characterise the responses in different tissue locations during collective electrotaxis. They found that the edges of the MDCK monolayer have a preferred direction of movement, which is outward of the tissue and is only observed within approximately 200 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} m of the tissue boundary. This effect is a hallmark of the collective migration of MDCK cells (Wolf et al. 2022; Tlili et al. 2018; Cohen et al. 2014; Alert and Trepat 2020). Although some attention has been given to understanding how the preferred migratory direction of the edges interferes with electrotaxis in the literature (Wolf et al. 2022), it is not necessary to incorporate detailed mechanisms of edge migration when understanding the average velocity of the centre of mass of the epithelial monolayer. Since the majority of the tissue does not exhibit the preferred direction of migration that the edges show, we restrict our analysis to the tissue bulk, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3\times 3\hbox { mm}^2$$\end{document} region in the centre of the tissue characterised by Wolf et al. (2022)—for a schematic, see Fig. 1. We take a spatial average of the velocity in the bulk, which yields one-dimensional data traces that represent the monolayer velocity in each of the experiments.
An Adaptation-Excitation Model of Collective Electrotaxis
We use a known model for cellular excitation and adaptation proposed by Erban and Othmer (2005) to model the internal cellular response to electric field. Since all the experiments we study consider a uniaxial electric field, we describe the migration of the monolayer bulk using the x component of the velocity, v, in the positive direction of the electric field, which is assumed aligned with the x-axis, as shown in Fig. 1. The goal of this section is to connect a model of cellular adaptation and excitation to the temporal evolution of the bulk velocity in response to the electric field and furthermore, to seek to use this model to address the question of how to optimally control the collective migration of MDCK epithelia by varying the field strength.
Model Derivation
We assume that cells experience a viscous friction force per unit mass which is proportional to their velocity with a friction coefficient, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} , such that the friction force per unit mass is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\gamma v.$$\end{document} The assumption of linear friction is well-founded for epithelia, since this corresponds to the polymerisation of adhesions on the substrate Alert and Trepat (2020). We assume throughout that there will be a uniaxial electric field in the positive v-direction whose magnitude is given by a scalar signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s = s(t)$$\end{document} . In what follows, we will introduce in detail the Erban and Othmer model Erban and Othmer (2005) to describe how the electric field creates an intracellular effective signal to the cytoskeletal machinery to produce an active force in the direction of the electric field.
The canonical model for excitation-adaptation proposed by Erban and Othmer considers an external stimulus, s(t), an effective signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}}$$\end{document} , and an inhibitor, I. The effective signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}}$$\end{document} , is a scalar quantity that describes the extent to which an external signal is transduced within the cell. There is therefore a difference between the constant electric field supplied and the resulting internal signal. The inhibitor, in contrast, represents the strength of the internal signalling pathways that try to inhibit this response to the external stimulus. The dynamics are given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \dot{s}_{\text {eff}}&= \frac{s(t) - (s_{\text {eff}} + I)}{\tau _e}, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \dot{I}&= \frac{s(t) - I}{\tau _a}. \end{aligned}$$\end{document}In this model, one assumes that the excitation timescale, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e > 0$$\end{document} , is far shorter than the adaptation timescale, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a > 0$$\end{document} . The different timescales in this system create a characteristic response to a constant signal Erban and Othmer (2005): when the external signal, s(t), is a Heaviside function, the solution of Eqs. (1) and (2) is of the form
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} s_{\text {eff}} = \frac{\tau _a}{\tau _a - \tau _e}\left( e^{-t/\tau _a} - e^{-t/\tau _e}\right) , \quad I = 1-e^{-t/\tau _a}. \end{aligned}$$\end{document}With other words, as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t \rightarrow \infty $$\end{document} , the effective signal wanes, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}} \rightarrow 0$$\end{document} . In particular, since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e \ll \tau _a$$\end{document} ,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} s_{\text {eff}} \sim \frac{\tau _a}{\tau _a - \tau _e} e^{-t/\tau _a}, \end{aligned}$$\end{document}whenever \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t \gg \tau _e$$\end{document} . This implies that on time scales larger than the excitation time, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} , the field response decays exponentially, so that a constant electric field will not produce an internal signal that leads to sustained migration at timescales much greater than the initial excitation time, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} . In this work, we consider s, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}}$$\end{document} , and I to be unitless, as we rescale the field strength, s, by dividing by the canonical field strength of 3V/cm in the experimental set-up. Finally, we make the assumption that whenever \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s \ge 0$$\end{document} , the effective signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}}$$\end{document} , equally satisfies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}} \ge 0$$\end{document} . The reason for making this additional assumption is that Eqs. (1) and (2) allow for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}}$$\end{document} to become negative when the field is turned off. MDCK epithelial monolayers are not observed to move in the direction opposite to the electric field when the electric field is turned off, and this is not physically realistic. For this reason, going forward, unless explicitly noted, we consider the quantity
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} s_{\text {eff}}^+ = \max (s_{\text {eff}}, 0). \end{aligned}$$\end{document}To close the system of equations, we now relate the system of effective signal and inhibitor to the velocity. We assume that the active force per unit mass in the direction of the electric field is proportional to the effective signal. This assumption follows canonical models in the literature for coupling polarity and force magnitude Alert and Trepat (2020), and results in the following expression for the active force per unit mass, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{F}_{\text {active}}$$\end{document} ,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \tilde{F}_{\text {active}} = \alpha s_{\text {eff}}^+, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha > 0$$\end{document} is a parameter that controls the responsiveness of the bulk to the input signal. Note that, since s is unitless, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} is given in units of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} m/ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {h}^2$$\end{document} . In this section, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} is assumed to be a constant, but we will explore different functional forms for the field response in the next section. By then assuming that the resultant force on the bulk is given by the sum of the active force per unit mass and the friction force per unit mass, one obtains a single equation for the velocity,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \dot{v} = -\gamma v + \alpha s_{\text {eff}}^+. \end{aligned}$$\end{document}We assume that the monolayer is stationary at the beginning of stimulation and the electric field is turned off, so that the initial condition is given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} v(0) = s_{\text {eff}}(0) = I(0) = 0. \end{aligned}$$\end{document}Linear Stability Analysis
We begin by noting that Eqs. (1), (2), and (7) can be written as a linear system of equations such that
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}}{\text {d}t}\begin{pmatrix} v \\ s_{\text {eff}} \\ I\end{pmatrix} = \begin{pmatrix} -\gamma &{} \alpha &{} 0 \\ 0 &{} -\tau _e^{-1} &{} - \tau _e^{-1} \\ 0 &{} 0 &{} -\tau _a^{-1} \end{pmatrix} \begin{pmatrix} v \\ s_{\text {eff}} \\ I\end{pmatrix} + \begin{pmatrix} 0 \\ \tau _e^{-1}s(t) \\ \tau _a^{-1} s(t)\end{pmatrix}, \end{aligned}$$\end{document}subject to the initial conditions in Eq. (8). Note that the eigenvalues of the matrix in the homogeneous part of the right-hand side of Eq. (9) are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\gamma , -\tau _e^{-1}, -\tau _a^{-1}$$\end{document} , making the matrix invertible. At steady state, one has that the left-hand side of Eq. (9) is zero, and that the steady state solution, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(v^\star , s_{\text {eff}}^\star , I^\star )$$\end{document} , for a constant field strength, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s(t) \equiv s$$\end{document} , is given by taking the inverse,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{pmatrix} v^\star \\ s_{\text {eff}}^\star \\ I^\star \end{pmatrix} = \begin{pmatrix} -\gamma ^{-1} &{} -\frac{\alpha \tau _e}{\gamma } &{} -\frac{\alpha \tau _a}{\gamma } \\ 0 &{} -\tau _e &{} \tau _a\\ 0 &{} 0 &{} -\tau _a \end{pmatrix} \begin{pmatrix} 0 \\ \tau _e^{-1} s \\ \tau _a^{-1} s \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ s \end{pmatrix}. \end{aligned}$$\end{document}Therefore, Eq. (9) has a single steady state. Since this steady state is given by zero velocity and zero effective signal, it can be interpreted as corresponding to a tissue that has fully adapted to the electric field and exhibits no response. Moreover, this steady state is stable since all eigenvalues of the matrix on the right-hand side of Eq. (9) are strictly negative. This explains the waning response with constant stimulation, as there is no steady state where a non-zero velocity is attained.
Bayesian Inference of Model Parameters
Our objective is to calibrate the simple adaptation-excitation model in Eq. (9) to the experimental bulk velocity data of Wolf et al. (2022) by performing Bayesian inference on the unknown model parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma , \alpha , \tau _e, \tau _a$$\end{document} from the available bulk data.
First, we note that the velocity decay parameter, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} , can be estimated directly from the data, since Eq. (7) predicts that, given an initial velocity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0$$\end{document} , the bulk velocity will decay exponentially in the absence of an electric field. Since the experimental data contain information on the bulk velocity after stimulation with the electric field is turned off, it follows that the observed velocity decay can be used to immediately infer the decay rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} . From these considerations, we seek to fit an exponential of the form
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} v(t) = C e^{-\gamma (t-t_{\text {end}})}, \end{aligned}$$\end{document}to the bulk velocity data in the 2.5 h after stimulation with the electric field is turned off. Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{\text {end}}$$\end{document} is the time at which the electric field is switched off, and C is given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} C = v(t_{\text {end}}). \end{aligned}$$\end{document}We perform a least-squares fit for the decay rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} , to give an estimate of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 1.765$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {h}^{-1}$$\end{document} . Figure 2 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} as fixed going forward. The benefit of directly estimating \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} from data is that it reduces the dimensionality of the parameter space, thus simplifying the inference of the remaining parameters.Fig. 2. 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 exponential decay in the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C e^{-\gamma t}$$\end{document} (solid blue line). The least-squares solution gives an excellent fit to the data and can be used to identify the decay rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} (color figure online)
We now keep \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 1.765\hbox {h}^{-1}$$\end{document} fixed and perform Bayesian inference of the remaining parameters, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} . We use Markov chain Monte Carlo (MCMC) with a Haario-Bardenet adaptive covariance and four chains under the assumption of a Gaussian error implemented in Python with the PINTS pacakge Clerx et al. (2019). For the distribution of the error, we assume Gaussianity since previous experiments by Wolf et al. demonstrate that the experimental variability around the mean is approximately Gaussian Wolf et al. (2022). The (standard) choice of four MCMC chains reflects the fact that each run of the MCMC algorithm is a random walk, so that the use of several random walks can be used to assess the ability of the algorithm to converge to the region in parameter space with high posterior probability Clerx et al. (2019). As a metric of convergence we use the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{R}$$\end{document} statistic Clerx et al. (2019), which is defined as the ratio of the between- and within-chain variances of the chains (Clerx et al. 2019). The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{R}$$\end{document} statistic summarises mixing and stationarity of the chains, with a value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{R} = 1$$\end{document} showing perfect mixing and stationarity of the chains. We use the reference value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{R} = 1.05$$\end{document} for four chains (Vehtarh et al. 2021) and set a maximum number of MCMC iterations to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\cdot 10^4$$\end{document} . We use as data the velocity profiles from the 3V/cm pulse experiments from Wolf et al. (2022). The MCMC procedure yields a well-identified posterior distribution, with marginal posterior means 149.92 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} m/ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {h}^2$$\end{document} , 0.260h, and 2.038h, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} , respectively. These posterior means are in good agreement with previous intuition that the time scale for adaptation is longer than the time scale for excitation, as well as with the data, suggesting that the bulk slows down on a timescale of the same order of magnitude as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} . Plots of the marginal posterior distributions of the model parameters are shown in Fig. 3. These show that the parameters can be confidently identified from the data, given their well-defined posterior distributions.
To assess the predictive capability of the model given the posterior distribution of the model parameters, we compute the posterior predictive intervals of the model, i.e. a 95% confidence interval for the model outputs given the posterior distribution of the model parameters. These posterior predictive intervals are shown in Fig. 3. It can be seen that the posterior predictive intervals are tightly concentrated around the mean, showing little model uncertainty in predicting the experimental data, given the posterior distribution of the model parameters. By overlaying the posterior predictive intervals for the model onto the experimental data, we observe that the model predictions are in good agreement with the experimental data and lie comfortably within the standard deviation of the experimental data. We conclude that the adaptation-excitation model presented in this section, when calibrated to the experimental data, is able to recapitulate the phenomenon of attenuating velocity during electrotaxis to the field when the electric field is held constant.Fig. 3. Bayesian inference on a simple model of adaptation and excitation given by Eq. (9) using experimental data from an electrotaxis experiment. Left: bulk velocity using a constant field of 3V/cm. Posterior mean of the model together with a confidence interval between the 5% and 95% quantiles of the model posterior predictions in cyan, experimental data and their confidence interval in grey. Notice the small width of the posterior predictive intervals. Right panels, from top to bottom: posterior distributions for model parameters, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e, \tau _a$$\end{document} , respectively (color figure online)
Optimal Control
Having established that excitation-adaptation dynamics can be used to describe the velocity of the tissue bulk during electrotaxis, we turn to the question of designing a stimulation protocol that is in some way optimal for the experiment. We previously established that the question of how to design electric field stimulation protocols has received little attention due to the difficulties in predicting tissue response a priori. Instead, experiments have been carefully calibrated through trial and error by the experimental laboratories that have performed them over the past few years. At this point, it is worth noting that these experiments are costly and time-consuming, making a grid search or other brute-force methods to determine optimal temporal stimulation patterns infeasible.
In this section, we will consider the problem of predicting an optimal pattern of stimulation using the adaptation-excitation model, drawing extensively on optimal control theory. For a legible introduction, we refer the reader to Lenhart and Workman (2007). We will consider that the optimal control problem can be cast as an optimization problem of a functional, E, that depends continuously on the velocity as well as the signal,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} E(v,s) = \int _0^T f(v,s) \text {d}t, \end{aligned}$$\end{document}and that there is a fixed start and end time of stimulation, given by 0, T, respectively. We need to penalise for the total amount of electric current administered during the experiment, as applying large currents throughout the cell medium for large amounts of time yields damage to the tissue and the experimental set-up. For now, we will incorporate this by requiring that a proxy of the total amount of charge delivered to the tissue during stimulation remains constant, i.e.
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _0^T s(t)^2 \text {d}t= s_{\text {budget}}. \end{aligned}$$\end{document}Later in this section, we will investigate the role of introducing lower and upper bounds on the field strength instead. Note also that in Eq. (14) we have introduced a quadratic to account for the fact that the electric field can be positively or negatively charged, and that fields of both signs will affect the tissue.
Under the condition of Eq. (14), we set up the first optimal control problem to maximise the total distance travelled by the tissue. This is a straightforward question to ask when one is concerned with increasing the controllability of collective migration. We therefore consider the optimization problem
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \max _{s} \int _0^T v(t)\text {d}t. \end{aligned}$$\end{document}This optimisation problem is then subject to the condition in Eq. (14). In optimal control theory, this is known as an isoperimetric problem. In this formulation, we introduce a new variable,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} z(t) = \int _0^t s(t)^2 \text {d}t, \end{aligned}$$\end{document}such that the problem can be described as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \max _s E&= \int _0^T v \text {d}t, \quad \text {subject to} \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \dot{z}(t)&= s(t)^2, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} v(0)&= 0, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} z(0)&= 0, \quad z(T) = s_{\text {budget}}. \end{aligned}$$\end{document}Recall that the isoperimetric contraint above has been introduced to enforce a physical limit on the amount of charge to which cells in the monolayer are exposed. One might also interpret the isoperimetric constraint as a manner of interrogating how, given a fixed amount of charge, one might optimally distribute the delivery of this charge so that the total distance travelled by the monolayer is optimal. To make direct comparisons with the experiments performed in Wolf et al. (2022), we wish to make the stimulation pattern have the same integrated charge for constant stimulation at 3 V/cm, so we set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {budget}} = 27$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {V}^2$$\end{document} h/ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^2$$\end{document} . For this problem the Hamiltonian is given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} H&= v + \lambda _v(t)\dot{v} + \lambda _{s_{\text {eff}}}(t) \dot{s}_{\text {eff}} + \lambda _I \dot{I} - \mu _z \dot{z}\\&= v + \lambda _v(t) \left[ -\gamma v + \alpha s_{\text {eff}}\right] \\&+ \tau _e^{-1}\lambda _{s_{\text {eff}}}(t)\left[ s - s_{\text {eff}} - I \right] + \tau _a^{-1}\lambda _I(t)\left[ s-I\right] - \mu _z s^2,\\ \end{aligned} \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial \lambda _\bullet /\partial t = - \partial H/\partial \bullet $$\end{document} , such that
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _v}{\text {d}t}&= - \frac{\partial H}{\partial v} = -1 + \gamma \lambda _v, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _{s_{\text {eff}}}}{\text {d}t}&= - \frac{\partial H}{\partial s_{\text {eff}}} = -\alpha \lambda _v + \tau _e^{-1} \lambda _{s_{\text {eff}}}, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _{I}}{\text {d}t}&= - \frac{\partial H}{\partial I} = \tau _e^{-1} \lambda _{s_{\text {eff}}} + \tau _a^{-1} \lambda _I, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\mu }{\text {d}t}&= -\frac{\partial H}{\partial z} = 0 , \end{aligned}$$\end{document}subject to the terminal conditions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _\bullet (T) = 0$$\end{document} . Additionally, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = \mu _\star $$\end{document} is a constant. The optimality condition is given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} 0 = \frac{\partial H}{\partial S} = \tau _e^{-1} \lambda _{s_{\text {eff}}} + \tau _a^{-1}\lambda _I - 2\mu _\star s. \end{aligned}$$\end{document}The solution to Eq. (26) yields the optimal stimulation pattern, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s^\star $$\end{document} , which is given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} s^\star = \frac{1}{2\mu _\star }(\tau _e^{-1} \lambda _{s_{\text {eff}}} + \tau _a^{-1}\lambda _I) = \frac{\tau _e^{-1}}{2\mu _\star }\left( \lambda _{s_{\text {eff}}} + \frac{\tau _e}{\tau _a}\lambda _I\right) = \frac{\tau _e^{-1}}{2\mu _\star }\left( \lambda _{s_{\text {eff}}} + \epsilon \lambda _I\right) , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon = \tau _e/\tau _a$$\end{document} . To find the optimal stimulation pattern, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s^\star $$\end{document} , one must seek to solve the linear inhomogeneous system of equations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _\bullet $$\end{document} . Note that Eq. (22), the equation for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _v$$\end{document} , decouples from the other adjoint equations, so we can solve directly to find
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda _{v}\left( t\right) = \frac{{\left( \gamma \lambda _{v}\left( 0\right) - 1\right) } e^{\gamma t}}{\gamma } + \frac{1}{\gamma }. \end{aligned}$$\end{document}Applying the terminal condition, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _v(T) = 0$$\end{document} , yields that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _v(0) = -(e^{-\gamma T}-1)/\gamma $$\end{document} , and so we have
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda _{v}\left( t\right) = (e^{\gamma (t-T)}-1)/\gamma . \end{aligned}$$\end{document}Direct substitution of the solution for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _v$$\end{document} into Eq. (23), the equation for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{s_{\text {eff}}}$$\end{document} , then yields
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \lambda _{s_{\text {eff}}}&= \frac{\alpha }{\gamma (\gamma - \tau _e^{-1})} - \frac{\alpha }{\tau _e^{-1}(\gamma - \tau _e^{-1})} + c_1 e^{\tau _e^{-1}t} - \frac{\alpha e^{\gamma (t-T)}}{\gamma (\gamma - \tau _e^{-1})}\\&= -\frac{\alpha }{\tau _e^{-1}\gamma } + c_1 e^{\tau _e^{-1}t} - \frac{\alpha e^{\gamma (t-T)}}{\gamma (\gamma - \tau _e^{-1})}. \end{aligned} \end{aligned}$$\end{document}We can now fit the terminal condition, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{s_{\text {eff}}}(T) = 0$$\end{document} , such that
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} 0 = \lambda _{s_{\text {eff}}}(T) = -\frac{\alpha }{\tau _e^{-1}\gamma } + c_1 e^{\tau _e^{-1}T} - \frac{\alpha }{\gamma (\gamma - \tau _e^{-1})}, \end{aligned}$$\end{document}which gives
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} c_1 = e^{-\tau _e^{-1}T}\left( \frac{\alpha }{\tau _e^{-1}\gamma } + \frac{\alpha }{\gamma (\gamma - \tau _e^{-1})}\right) = \frac{\alpha \gamma }{\tau _e^{-1} (\gamma - \tau _e^{-1})}e^{-\tau _e^{-1}T}. \end{aligned}$$\end{document}Therefore,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda _{s_{\text {eff}}} = \frac{\alpha \gamma }{\tau _e^{-1}(\gamma - \tau _e^{-1})} \left( e^{\tau _e^{-1}(t-T)} - \frac{\tau _e^{-1}}{\gamma ^2}e^{\gamma (t-T)} \right) -\frac{\alpha }{\tau _e^{-1}\gamma }. \end{aligned}$$\end{document}We can then find \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _I$$\end{document} by solving Eq. (24),
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda _I = e^{-\tau _a^{-1}t}\int _0^t \lambda _{s_{\text {eff}}}(\xi ) e^{\tau _a^{-1}\xi }\text {d}\xi + C_I e^{\tau _a^{-1}t}, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_I$$\end{document} is determined by setting the terminal condition, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _I(T) = 0$$\end{document} , which gives
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} C_I = -e^{-\tau _a^{-1}T}\int _0^T \lambda _{s_{\text {eff}}}(\xi ) e^{\tau _a^{-1}\xi }\text {d}\xi . \end{aligned}$$\end{document}Therefore,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} s^\star = \frac{\tau _e^{-1}}{2\mu _\star } \lambda _{s_{\text {eff}}} + \epsilon \frac{\tau _e^{-1}}{2\mu _\star }\left[ e^{-\tau _a^{-1}t}\int _0^t \lambda _{s_{\text {eff}}}(\xi ) e^{\tau _a^{-1}\xi }\text {d}\xi + C_I e^{\tau _a^{-1}t} \right] . \end{aligned}$$\end{document}To determine the solution of the system, it remains to find the constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _\star $$\end{document} such that the optimal solution integrates to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {budget}}$$\end{document} . We use the marginal posterior means for the model parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha , \tau _e, \tau _a$$\end{document} to first compute \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s^\star $$\end{document} numerically and use numerical integration to determine the integral during the stimulation window.Fig. 4. Stimulation patterns for optimal control of collective electrotaxis. In all plots, red corresponds to maximum tissue displacement (optimal solution given by Eq. (27)), and yellow corresponds to maximum terminal velocity (optimal solution given by Eq. (32)). Blue curves represent 3V/cm electric field pulse stimulation pattern used experimentally. Left: the different stimulation patterns in normalized units. Middle: corresponding velocity curves using the different stimulation patterns. Right: data from 3V/cm pulse experiment, compared to posterior distribution for optimal solution for distrance travelled and terminal velocity. Parameter values correspond to the posterior means for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} , and the pre-computed value for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} (color figure online)
We show a prediction for the optimal stimulation pattern in Fig. 4. Interestingly, the optimal stimulation pattern begins at a nonzero electric field value and increases for the first half of the stimulation window, when it exceeds the constant electric field strength of 3V/cm, before dropping below that field strength in the second half of the stimulation window. The optimal stimulation pattern found in Eq. (27) shows only a modest increase in the distance travelled by the bulk of the monolayer compared to the experimental data, since the integrated numerical solution of the velocity is 2.17% greater than the distance travelled using a constant electric field. This reflects the fact that electrotaxis experiments have been iterated many times experimentally. Strikingly, the optimal stimulation pattern also results in a bulk velocity that still exhibits a slowdown during stimulation. Commonly, studies on collective migration induced by electrotaxis aim to optimise electrotaxis by maximising the velocity at which cells migrate (Wolf et al. 2022; Leal et al. 2023; Cohen et al. 2014). If maximising the speed of collective migration is the goal, optimising collective migration during electrotaxis would be incompatible with allowing a significant slowdown during collective migration. Our findings, however, suggest that the slowdown of the bulk during stimulation does not immediately represent nonoptimality of the stimulation procedure, if maximising the total distance travelled instead is the aim.
We incorporate parameter and prediction uncertainty by computing the posterior distribution over the optimal control solution. We take \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\cdot 10^3$$\end{document} samples from the MCMC posterior distribution and for each sample compute the optimal stimulation protocol according to Eq. (27) using that given parameter combination and recomputing the normalising constant, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _\star $$\end{document} , for each new parameter combination to obtain a normalised density. We again compute the 5% and 95% percentiles of the velocity distribution at each time point and plot them in Fig. 4. We see that the posterior distribution is tight around the posterior mean, which confirms that the above analysis holds true across the entire posterior distribution.
Optimal Control for the Terminal Velocity
Optimising for the total distance travelled predicts that there is only a very moderate improvement in the optimised outcome compared to stimulation with a constant electric field of 3V/cm. A striking feature of the optimal stimulation curve shown in Fig. 4 is that the velocity curve under the optimal stimulation protocol, like the pulse stimulation protocol, displays a diminishing response to the field at later times. Therefore, we pose the additional question of how to maximise the terminal velocity of the tissue, in other words: How can the stimulation be designed so that the bulk velocity at the terminal time, T, is maximal? A priori, one would expect this to be a different stimulation pattern than the optima found previously. We change the objective function so that we maximise the terminal velocity, v(T),
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} E(v) = \int _0^T \dot{v} \text {d}t = \int _0^T [-\gamma v + \alpha s_{\text {eff}}] \text {d}t. \end{aligned}$$\end{document}Note that this change in objective function only changes the equations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {d}\lambda _v/\text {d}t$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {d}\lambda _{s_{\text {eff}}}/\text {d}t$$\end{document} in the adjoint equations, Equations (22)-(25), since the only additional terms introduced into the Hamiltonian include the velocity, v, and the effective signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}}$$\end{document} . We obtain
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _v}{\text {d}t}&= \gamma + \gamma \lambda _v, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _{s_{\text {eff}}}}{\text {d}t}&= -\alpha \lambda _v + \tau _e^{-1}\lambda _{s_{\text {eff}}}, \end{aligned}$$\end{document}such that, directly fitting the terminal condition, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _v(T) = 0$$\end{document} , we have
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda _v = e^{\alpha (t-T)} -1, \end{aligned}$$\end{document}which we substitute into Eq. (23) to find
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda _{s_{\text {eff}}} = \frac{\alpha }{\alpha -\tau _e^{-1}}\left( e^{\tau _e^{-1}(t-T)} - e^{\alpha (t-T)}\right) . \end{aligned}$$\end{document}The same computation as in Eq. (27) now yields predictions of the optimal electric field stimulation protocol, which results in a terminal velocity of 62.9 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} m/h compared to 52.5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} m/h for a constant field. The electric field that optimises the terminal velocity is shown in Fig. 4. Note that the field strength is monotonically increasing to its maximum intensity at the end of stimulation. We obtain a higher terminal velocity than in the other experiments, at the cost of a very high field strength at the end of stimulation, approximately 6 V/cm, which is physiologically and experimentally feasible, albeit not for long periods of time. The fact that the optimal control problem for the terminal velocity yields a ramping pattern suggests that total distance travelled and terminal velocity are two incompatible requirements to fulfill simultaneously, given the model and evidence presented. Finally, to address parameter and prediction uncertainty, we repeat the computation of posterior distributions over the velocities, as done in Sect. 4, which we display in Fig. 4. Again, the posterior distribution of the velocities is tight around the posterior mean.
Optimal Control for Constant Velocity
The optimal control problems for the maximum distance travelled and the maximum terminal velocity seem to suggest that maximising the distance comes at the cost of a decrease in velocity during stimulation, while optimising the terminal velocity does not maximise the total distance travelled. The finding that the optimal stimulation pattern for total distance travelled does not yield a nondecreasing bulk velocity during the course of stimulation with the electric field runs counter to intuition in the electrotaxis community, as the current opinion in the field tends to view a drop in bulk velocity as a source of nonoptimality in stimulation. From an experimental viewpoint, controlling field strength so that the migratory velocity is constant is attractive, since electric fields are often used for ‘cruise control’. Therefore, we set up a final optimal control problem designed to keep variations in the velocity minimal. We do this by once more modifying the functional so that we seek solutions that minimise
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} E(\dot{v}) = \int _0^T \dot{v}^2 \text {d}t = \int _0^T (-\gamma v + \alpha s_{\text {eff}})^2\text {d}t . \end{aligned}$$\end{document}As in the previous section, the only adjoint equations affected are those for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _v$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{s_{\text {eff}}}$$\end{document} . These become
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _v}{\text {d}t}&= 2\gamma (- \gamma v + \alpha s_{\text {eff}}) + \gamma \lambda _v, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _{s_{\text {eff}}}}{\text {d}t}&= -2\alpha (-\gamma v + \alpha s_{\text {eff}}) - \alpha \lambda _v + \tau _e^{-1} \lambda _{s_{\text {eff}}}. \end{aligned}$$\end{document}In addition, we impose a terminal condition for the velocity, so that we prescribe a desired velocity at the end of stimulation, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v(T) = v^\star $$\end{document} for some pre-determined constant, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} , which can be set by the practitioner in accordance with a practical need. Together with Eqs. (9), (34), (35) and the usual transversality conditions, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _\bullet (T) = 0$$\end{document} , the optimisation of the function in Eq. (33) creates a nonlinear first-order boundary value problem (BVP). The complicated coupling between the equations requires a numerical solution to the problem. Here, we set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} so that it is equal to the average bulk velocity during stimulation with a single pulse of 3 V/cm.
For numerical solution of the BVP, we use the bvp4c scheme of Kierzenka et al. implemented in Python through the SciPy package Virtanen et al. (2020). This scheme implements a fourth-order collocation algorithm with control of the residuals and uses a damped Newton method with an affine-invariant criterion function. For the implementation, we set a maximum relative tolerance of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-3}$$\end{document} , an absolute boundary value tolerance of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-6}$$\end{document} and specify a maximum of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} nodes for the collocation algorithm.Fig. 5. Solution to optimal control problem for constant bulk speed (optimal solution for the problem in Eq. (33)) during stimulation with an electric field. Left: bulk speed for constant velocity (solid line) plotted alongside experimental data with constant electric field stimulation protocol (dashed line). Middle: numerical solution of the optimal control problem for constant velocity (solid line) plotted against reference normalized electric field strength (dashed line). Right: effective signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}}$$\end{document} . Parameter values correspond to the posterior means for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} , and the pre-computed value for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}
Figure 5 shows that the electric field profile needed to produce a constant velocity profile is very strong initially and rapidly decays to create a constant effective signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\text {eff}}$$\end{document} . In fact, the field strength predicted here is 75 V/cm, which is not physiologically realistic, as it would lead to toxicity and rapid cell death. This finding suggests that given the physiological constraints of the problem, it is not possible to design an optimal control experiment that creates a constant speed profile with the same average velocity as that arising from pulse stimulation. This large initial peak in electric field strength arises from the requirement that the tissue hits a constant velocity quasi-instantly upon stimulation. This finding invites the question of which constant velocity values are feasible to achieve given that the solution profile needs to be realistic, i.e. given the physical constraints on the electric field strength. To this end, we perform a grid search, whereby we solve the optimal control problem for a range of target velocities, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} , ranging from 0 to 40 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} m/h. For each of these terminal velocities, we solve the BVP for the optimal stimulation protocol and compute the resulting bulk velocities and optimal stimulation strategy. We summarise the results in Fig. 6.Fig. 6. Stimulation profiles and electric field strengths for varying terminal velocities in optimal solutions to the problem in Eq. (33)). From left to right. Bulk velocity when using optimal control solution for constant stimulation using different values of terminal velocity (see legend). Optimal stimulation profile corresponding to the different terminal velocities in the left panel. Electric field strength at the start of stimulation. Electric field strength averaged during stimulation. Note the difference in the magnitude of the axes. Parameter values correspond to the posterior means for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} , and the pre-computed value for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} (color figure online)
Figure 6 shows that at each of the values for the terminal velocity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} , the model predicts that an optimal stimulation strategy can be designed so that the bulk moves at this constant speed. Each of these stimulation profiles has the same shape, whereby a very large current is applied at the beginning of the experiment, which slowly decreases during the course of the experiment. We investigate the relationship between the desired terminal velocity and the field strength at the onset of stimulation as required by the model. We see that there is a linear increase in the normalised electric field density with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^{\star }$$\end{document} and that the initial electric field strength is prohibitively large for velocities higher than 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document} m/h. This is surprising and experimentally frustrating given that the mean electric field strength predicted by the model only shows a marginal increase with the desired terminal velocity, whereby all terminal velocities experience a normalised electric field strength that is physiologically tolerable. We conclude, thus, that designing a realistic electric field to maintain the tissue bulk at constant velocity, while keeping the electric field strength within realistic bounds, is not possible.
The source of the unphysically strong electric field at the start of stimulation is the requirement that the field immediately reaches the desired velocity. At the same time, the velocity in the period after the onset of stimulation is constant, with a far smaller field strength than that required at the beginning. Motivated by these two facts, we ask if we can overcome the problem of non-realistic predictions for the initial field strength by setting a maximum field strength allowed for s. We expect this might come at the cost of reaching the target velocity at the beginning of stimulation. We allow for the field to have a field strength of at most 9V/cm, which is around the limit at which electric fields are tolerated in vitro. For the optimal control problem in Eq. (33), we require that the physical optimal solution, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s^\star _{\text {phys}}$$\end{document} , satisfies
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} s_{\text {phys}}^\star = \min \left( s^\star , s_{\max }\right) , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\max }$$\end{document} is a parameter we can vary. By considering \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\max }$$\end{document} across a range between 3 and 9 V/cm, we solve the optimal control problem and show the results in Fig. 7. We set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} as the mean velocity in the 3V/cm pulse experiment and now vary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\max }$$\end{document} . While the solutions are qualitatively similar and the velocity is constant during most the experiment, the solutions are not identical. A closer inspection of the terminal velocity error shows that solutions with a weaker initial field have a larger error at the end, implying that the large initial pulse is what is necessary to create the terminal velocity, but it is the temporal evolution of the field which is necessary to create a plateau. While this error decreases rapidly for small values of the maximum field strength, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{\max }$$\end{document} , a close agreement with the desired cruise velocity is only achieved, once again, with unfeasibly high electric field strengths.Fig. 7. Optimal stimulation to keep the velocity constant during stimulation. Field strengths are solutions are solutions to the problem in Eq. (33) with maximum field strength defined as in Eq. (36). Left: bulk velocity using different maximum electric field strengths (see legend). Middle: physical electric field profile showing similar stimulation protocols. Right: terminal velocity error decreases when the maximum normalized electric field strength increases. Parameter values correspond to the posterior means for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} , and the pre-computed value for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} (color figure online)
Constant Velocity During a Small Time Interval
When the goal was to optimise to keep the velocity constant, as in the previous section, we observed that the electric field profile had a non-physically high field density at the start of stimulation, making this optimal control strategy infeasible in practise. This comes from requiring that the velocity hits the target velocity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} , quasi-instantly upon the start of stimulation. Here, we ask the question of how solutions change when instead the velocity is only required to be constant on an interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[t_1, t_2]$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0< t_1< t_2 < T$$\end{document} . Here, we will explore how the optimal electric field profile is affected by this requirement. We will also study how solutions are affected by altering the values for the start and end of stimulation, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_2$$\end{document} , respectively.
The functional in Eq. (33) needs to be modified to address the optimality condition that the velocity is constant on the interval. This is because minimising the square of acceleration \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{v}^2$$\end{document} as in Eq. (33) leads to a trivial solution for the velocity: when the velocity is assumed to be zero by the initial conditions and the velocity is continuous, the optimal solution will be zero throughout the stimulation time frame. Therefore, one must consider a functional that directly relates the velocity during the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[t_1, t_2]$$\end{document} to the target velocity. Such a functional can be defined as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} I(v) = \int _0^T \mathbb {I}(t \in [t_1, t_2])\cdot (v-v^\star )^2 \text {d}t, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {I}$$\end{document} is an indicator function. This indeed restricts the constant velocity requirement to the desired interval, but the integrand is no longer a continuous function of time. To make the optimisation problem admissible, we choose to approximate the indicator function in the integrand by a smooth bump function, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi $$\end{document} , which we define as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Phi (t) = \frac{1}{1 + \exp \left( - \frac{t - t_1}{\epsilon }\right) }\cdot \frac{1}{1 + \exp \left( - \frac{t_2 - t}{\epsilon }\right) }, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} determines the length scale of the region where the function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi $$\end{document} transitions from zero to one. We set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon = 50^{-1}$$\end{document} so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \ll t_1$$\end{document} in all of our considered optimisation problems. Now we obtain an admissible optimal control formulation, given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} E(v) = \int _0^T \Phi (t)\cdot (v-v^\star )^2 \text {d}t, \end{aligned}$$\end{document}which yields a modified version of the adjoint equations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _v$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{s_{\text {eff}}}$$\end{document} , which are Eqs. (34) and (35), respectively. These become
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _v}{\text {d}t}&=\Phi (t)\cdot 2(v - v^\star ) + \gamma \lambda _v, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\text {d}\lambda _{s_{\text {eff}}}}{\text {d}t}&= - \alpha \lambda _v + \tau _e^{-1} \lambda _{s_{\text {eff}}}, \end{aligned}$$\end{document}subject to the terminal conditions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _\bullet (T) = 0$$\end{document} . Note that, in contrast to the optimal control problem for constant velocity in the previous subsection, there is no need to set up a terminal condition for the velocity. This means that the zero initial conditions of Eq. (8) can be imposed. As before, we solve this boundary value problem using the bvp4c numerical scheme introduced in Sect. 4.2.
To determine the range of admissible values for the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[t_1, t_2]$$\end{document} , we make the assumption that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_2 = T-t_1$$\end{document} , i.e. the stimulation interval is symmetric within the interval [0, T]. A priori, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} should be chosen sufficiently large that a physically realistic electric field is optimal, but not so large that the field is off at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0$$\end{document} . To understand how to choose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} , then, we note that there exist two characteristic times in this optimal control problem. First, since constant stimulation for times longer than the adaptation time scale, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} , results in a waning response from the field, choices for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} that result in the stimulation time being too long, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(T-t_1)/2 \gg \tau _a$$\end{document} , will require a large value of the electric field initially to prevent cellular adaption to the field and will hence provide unphysical solutions that are not realistic stimulation protocols in practice. At the same time, one might expect that there is a time scale that dictates when the field is initially on at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} h and when it is off: if the waiting time from the beginning of stimulation, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} , is shorter than this critical time scale, the electric field should be turned on immediately, since otherwise the monolayer will not reach the target velocity during the optimality time window. Intuitively, this time scale should be of the same order as the time it takes for the velocity to reach its maximum under a constant electric field. We can estimate of this time scale, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\max }$$\end{document} , by maximising the expression in Eq. (3), which yields
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \tau _{\max } = \frac{\log (\tau _a) - \log (\tau _e)}{\tau _e^{-1} - \tau _a^{-1}}. \end{aligned}$$\end{document}For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1 > \tau _{\max }$$\end{document} one would therefore expect that the optimal electric field density is zero at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0$$\end{document} , while for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1 < \tau _{\max }$$\end{document} one would expect that the density is nonzero. Therefore, one would expect that choices of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} in the region
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} t_1 \in \left[ \frac{T-\tau _a}{2}, \tau _{\max }\right] , \end{aligned}$$\end{document}will yield optimal electric field protocols that are high at the onset of stimulation when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} is at the left boundary of the interval, or electric field protocols that slowly ramp up after the onset of stimulation. We also expect values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} smaller than the values in this range to yield physically unrealistic electric field protocols, and solutions for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} larger than the values in this range to have a near-zero initial electric field strength. In Fig. 8, we show how the optimal electric field protocol changes when different values for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1, t_2, v^\star $$\end{document} are used.Fig. 8. Optimal control to keep velocity constant in a small time window (solutions to the optimal control problem in Eq. (39)). Left: bulk velocity traces using different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} for stimulation. Middle: electric field strengths with different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} . Right: initial electric field strength in the optimal control solution as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} . Parameter values correspond to the posterior means for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _e$$\end{document} , and the pre-computed value for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} (color figure online)
Figure 8 shows that the optimal electric field protocols change qualitatively as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} is varied. For small values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} , close to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} , the optimal electric field protocol exhibits a sharp overshoot of the velocity prior to reaching a plateau in the desired interval. This is achieved by the optimal electric field protocols exhibiting a sharp peak initially, bringing the velocity above the target, and then letting the velocity decrease to the target. This overshoot disappears whenever \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1 > \tau _{\max }$$\end{document} . In this regime, the field allows the velocity to slowly increase to the desired target velocity. All optimal electric field protocols exhibit a ramp. This effect is clearly visible in the final plot of Fig. 8, which shows that there is a sharp decrease in the initial field strength whenever \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1 < \tau _{\max }$$\end{document} . We conclude that a constant velocity can be achieved without having an overshoot when the waiting time, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1$$\end{document} , is sufficiently long. Finally, we note that by performing a grid search across possible values for the target velocity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} , we found that these results do not change qualitatively with the value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} , as the only change is that the maximum magnitude of the optimal electric field protocol scales linearly with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^\star $$\end{document} .
Comparison to Bang–Bang Control
The modest gain in distance travelled by the tissue predicted when the optimal stimulation protocol is used suggests that naive stimulation by keeping the electric field at constant field density might be an effective policy. This is not surprising given that it is the policy converged upon after years of experimental iteration, and that the key variable here is when to switch off the field. This is reminiscent of bang-bang control in linear optimisation theory. Bang-bang control occurs when the Hamiltonian, H, of the optimal control problem is linear in the control variable, so that it disappears from the optimality condition altogether (Lenhart and Workman 2007). In our problem the Hamiltonian is nonlinear due to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _zS^2$$\end{document} term included to enforce the condition that a fixed amount of charge is injected into the system. If one now assumes instead that the field is always positive, the condition could be relaxed so that instead S is bounded by some maximum strength, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 \le S \le S_{\max }$$\end{document} . When the isoperimetric constraint is removed, the Hamiltonian becomes
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} H = v^2 + \lambda _v(t) \left[ -\gamma v + \alpha s_{\text {eff}}\right] + \tau _e^{-1}\lambda _{s_{\text {eff}}}(t)\left[ s - s_{\text {eff}} - I \right] + \tau _a^{-1}\lambda _I(t)\left[ s-I\right] ,\nonumber \\ \end{aligned}$$\end{document}with optimality condition
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial H}{\partial S} = \tau _e^{-1} \lambda _{s_{\text {eff}}} + \tau _a^{-1}\lambda _I, \end{aligned}$$\end{document}so that the optimal control is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s^\star = s_{\max }$$\end{document} , whenever the right-hand side of Eq. (45) is positive by the Pontryagin optimality criterion. This, we have established in the previous numerical experiments, holds whenever \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t < T$$\end{document} . In sum, we find that naive stimulation with a pulse is the optimal stimulation policy when the condition of total injected charge is relaxed and replaced instead with an upper bound on the electric field strength.
Discussion and Outlook
In this work we have developed and calibrated a canonical model for excitation-adaptation to describe the cellular response to collective electrotaxis, and we have identified that the mechanism of excitation and adaption can faithfully describe the temporal dynamics of the bulk response to electric field stimulation. We have used this model to optimise for different experimental outcomes, including maximising the total distance travelled by the tissue bulk, maximising the terminal velocity, and keeping the migration velocity constant.
Our model shows that maximising the distance travelled by the tissue bulk is possible by changing the temporal distribution of electric field strength, but it only produces a marginal increase in performance relative to stimulation with an electric field with constant field strength. In fact, this increase is well within the standard deviation of the experimental data. We related the near-optimal performance of the constant stimulation experiment to the fact that constant stimulation is the solution to a bang-bang control problem. A similar insight arises from controlling the maximum terminal velocity of the tissue, which recovers a ramp, but has an unphysically large field strength at the end of electric field stimulation. Together, this challenges the idea in the electrotaxis community that a bulk velocity that is not constant during stimulation – but rather exhibits a sharp peak and a rapid decrease during stimulation –is not optimal, as this velocity profile is in close agreement with the optimal solution for the maximum distance travelled for the tissue bulk.
Most importantly, our analysis shows that it is not possible to generate a constant velocity profile during electrotaxis when the experimental time exceeds the adaptation time scale, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _a$$\end{document} . This is reflected in the prediction that an unphysical field strength is necessary to create a constant velocity profile for the duration of the experiment. When we relaxed the constraint to allow the tissue to slowly increase to the desired velocity and then remain constant, such optimal electric field protocols only yielded velocity profiles without an overshoot whenever the stimulation time satisfied certain constraints. Together, these findings show that there exist physical limits to how we can stimulate, both in terms of duration of the experiment, and physiological constraints on electric field strength.
In the context of collective electrotaxis, this work can be extended and applied in several possible ways. Firstly, our framework establishes a reproducible and predictive model for collective cell migration that can be used to develop new experimental approaches. Conversely, experimental data can help validate and develop the optimal control framework proposed in this work. More generally, our modelling approach can be generalised and extended to be used in other contexts where there exists a desired collective migration outcome and one must find an optimal electric field protocol to achieve this pattern of migration. This can be achieved directly in our framework by expressing the desired movement pattern as the optimal solution to a linear control problem and using the optimal control framework put forth in this work. Secondly, in this work we have only considered uniaxial electric fields, while experimental set-ups exist to apply spatially dependent fields (Cohen et al. 2014; Leal et al. 2023; Song et al. 2021). Our framework can be extended to optimise for desired migration outcomes that vary in both time and space. For the optimisation and control of spatially varying fields, one must first develop a spatial model of collective electrotaxis that takes into account the different migratory cues that are present in the different parts of the monolayer. Such a spatially resolved model will enable the design and implementation of physically relevant applications of electrotaxis, such as those used in wound healing or steering collective migration. Finally, this work provides a blueprint for using optimal control theory in conjunction with external inputs that can be modulated in time, to control cellular collectives. Such a tunable external input might, for example, be implemented in the context of chemotaxis, which was the external signal originally modelled by Erban and Othmer (2005). By carefully setting up an optimal control problem corresponding to the maintenance and adaptation of a chemical gradient, one could aid the experimental design of microfluidics or other ways of maintaining chemical gradients to ensure desired outcomes during chemotaxis.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Alert R Trepat X Physical models of collective cell migration Annu Rev Condens Matter Phys 20201117710110.1146/annurev-conmatphys-031218-013516 · doi ↗
- 2Allen GM Mogilner A Theriot JA Electrophoresis of cellular membrane components creates the directional cue guiding keratocyte galvanotaxis Curr Biol 201323756056810.1016/j.cub.2013.02.04723541731 PMC 3718648 · doi ↗ · pubmed ↗
- 3Clerx M Robinson M Lambert B Lei CL Ghosh S Mirams GR Gavaghan DJ Probabilistic inference on noisy time series (PINTS)J Open Res Softw 201971610.5334/jors.252 · doi ↗
- 4Cohen DJ Nelson WJ Maharbiz MM Galvanotactic control of collective cell migration in epithelial monolayers Nat Mater 201413440941710.1038/nmat 389124608142 · doi ↗ · pubmed ↗
- 5Erban R Othmer HG From signal transduction to spatial pattern formation in. E coli: a paradigm for multiscale modeling in biology Multisc Model Simul 20053236239410.1137/040603565 · doi ↗
- 6Gao R Zhao S Jiang X Sun Y Zhao S Gao J Borleis J Willard S Tang M Cai H Kamimura Y Huang Y Jiang J Huang Z Mogilner A Pan T Devreotes PN Zhao MA large-scale screen reveals genes that mediate electrotaxis in Dictyostelium discoideum Sci Signal 201510.1126/scisignal.aab 0562 PMC 447047926012633 · doi ↗ · pubmed ↗
- 7Leal J Shaner S Jedrusik N Savelyeva A Asplund M Electrotaxis evokes directional separation of co-cultured keratinocytes and fibroblasts Sci Rep 20231311144410.1038/s 41598-023-38664-y 37454232 PMC 10349865 · doi ↗ · pubmed ↗
- 8Li L He Y Zhao M Jiang J Collective cell migration: implications for wound healing and cancer invasion Burns Trauma 201311212610.4103/2321-3868.11333127574618 PMC 4994501 · doi ↗ · pubmed ↗