Multicellular adaptation to electrophysiological perturbations analyzed by deterministic and stochastic bioelectrical models

When we attempt to simulate the cell bioelectrical adaptation to depolarization caused by a potassium channel blockade, the broad repertoire of channels available makes the inverse problem difficult for cells: which of the different conductance combinations would counteract the externally-induced depolarization? However, the time required to explore every channel combination would be exceedingly large. Also, the number of possible solutions is severely restricted because channel and gap junction proteins have other physiological functions in addition to regulating the cell potential5,16,26–29,44. Thus, most of the potential combinations would not meet all homeostatic requirements and could then be dead-ends for the cell. In the above context, we have assumed that the membrane potential is the master regulation of cell bioelectricity and considered the forward problem of restoring its correct value using an effective res channel. To this end, we have introduced a voltage-gated feedback for the adaptation that eventually upregulates the expression of the res channel, which converts the inverse problem of giving precise values to multiple channel conductances into the forward problem of controlling the physiological cell potential with a reduced number of channels.

We must admit that many error signals and feedbacks, together with other master regulators, might be involved in real cases, but we could expect that the number of cell properties to be controlled should be lower than the total number of available channels. Thus, it can be advantageous to replace the inverse problem of establishing precise individual levels of channel expression by the hopefully simpler problem of identifying the small number of master regulators that control the relevant homeostatic properties. In this scenario, the membrane potential, which can be regulated by a limited number of counteracting channels16,19,20,33, could be a choice because of its central role in metabolic, osmotic, developmental, and regenerative processes5,19,20,27,31,34.

Here, the following charateristics of the bioelectrical cell potential5,16,27,34,45–47 should be emphasized: (1) it shows a high sensitivity to a broad range of chemical, mechanical, and electrical signals, responding rapidly to external changes in the environmemt; (2) it can be tuned by a variety of ions and charged molecules inside and outside the cell; and (3) it allows a multicellular scaling in terms of electric potential patterns that are morphologically instructive. Note also that the membrane potential provides an experimental control of the correlations among channel mRNA levels not only in non-excitable cells5,36 but also in neurons9,38. In this particular context, the channel and junction proteins may act as the cellular hardware while the cell potential could be the software layer that fixes the hardware layer problem –the blocking of the potassium channel5,16. In this way, the bioelectric software enables flexibility to the multicellular aggregate, allowing to respond to environmental stimuli without having to change the hardware.

To characterize the deterministic results further, we have assumed that the depolarization stressor makes the cell to relax the initially precise control of gene expression, thus allowing stochastic changes in the transcription of the res channel that permits the adaptation. The transition from tight control to plasticity is a universal adaptation mechanism characteristic of complex systems, from protein structural states to social network responses. The exploration algorithm used leads to results in qualitative agreement with the deterministic model, thus suggesting that the general concepts introduced can be useful.

Experimentally, the simulation results can be relevant to model animals such as planaria. Voltage-sensitive dyes reveal that distinct anterior-posterior planaria morphologies can be obtained by changing the axial electric potential24,25, which is regulated by the cell ion channels. In particular, the head degeneration of planaria (Dugesia japonica) observed after barium chloride exposure4 is ascribed to the Ba2+-induced blocking of potassium channels, which provokes the head depolarization4. Following this local depolarization, however, long-term changes in the channels expression of the head cells can lead to repolarization and counteraction of the planaria excitotoxic stress. Also, prolonged Ba2+ contact eventually results in the regeneration of a new head which maintains form and function despite the presence of barium. Since the neural networks of the brain are mostly destroyed after the first barium chloride exposure, they cannot be the only responsible factor for the adaptation, although they may still influence the subsequent regeneration via additional feedback loops.

In an electrophysiologically-focused description of the above adaptation process4, TRPM channels and possibly other alternative cation channels could be upregulated in BaCl2-exposed animals: the subsequent blocking of these channels prevents planaria adaptation4, as suggested by the model of Fig. 1. Note also that the system response depends on the presence of barium: when a normal water environment is re-established, the unblocking of the potassium channel resets the cell potential to the normal polarized value, thus suppresing the cell stress and then decreasing the rescue channel transcription. As a consequence, it is conceivable that the barium resistance decreased after a sufficiently long period in normal water, as observed experimentally4. The above facts suggest that the planaria adaptation could be partly attributed to an electrophysiological exploratory dynamics.

Using deterministic and stochastic models, we have ascribed the cell adaptation to the short-term bioelectrical perturbation that triggers long-term downstream biochemical processes because of the local correlation between the electric potential and the distribution of signaling ions and molecules that influence transcription. The basic concepts invoked here provide a complementary view to traditional biochemical models and give new insights into sensing and adaptation mechanisms that do not rely on neural networks. In the biophysical model, the cell does not need to fine-tuning the transcription of every channel protein in order to fix precise conductance values. Instead, it can make use of the cell potential5,16 in an error-correcting mechanism that involves a small subset of the total number of cell channels. In principle, these subset of channels can be experimentally identified by external silencing and blocking methods4,5. Note also that exploratory adaptations to environmental changes based on modularity48 and collectivity49 concepts that do not rely exclusively on fine-tuning molecular processes should make easier corrective actions in real systems50. Understanding the mechanisms of such error corrections is critical not only to the biomedicine of electrophysiological states, but also to the evolutionary biology of plasticity in development and regeneration as these processes adaptively handle highly novel scenarios51,52. In particular, local electric circuits based on ion channels and gap junction networks are crucial in collective cell migration and tissue morphogenesis processes44,53,54. Also, ion channels and cell potentials enable active communication and collective responses in bacterial communities55.

Here, we implicitly assume that responses to a general class of stressors (cell depolarizing agents here) can be elicited by a novel particular member of that class: while barium might be new to planaria4, its effect recalls previous depolarization events that should be counteracted. In this context, the superfamily of voltage-gated ion channels can offer an evolutionary ancient response mechanism17. At the single-cell level, it is plausible that the adaptive re-polarization be evolutionary encoded by the upregulated transcription of a rescue channel. At the multicellular level, it is suggested here that the adaptation can also involve a population-level process where the community effect of the polarized cells in the right region of the aggregate acts to counteract the depolarization process in the left region, thus restoring the bioelectrical pattern. Such non-cell-autonomous responses have been seen before in examples of defect repair56 and long-range bioelectric signaling in tumor formation57, and offer tantalizing targets as biomedical targets.

Figure 5 shows the equations of the biophysical model used in the simulations. Note that we do not propose simple equivalent circuits here but explore instead a complete description for the biological interplay between the genetic and the bioelectric networks. This description is extended further from the single-cell to the multicellular scale by the intercellular gap junctions16.

Initially, the cells have the isolated membrane potential value and then the multicellular system is allowed to relax to the physiological potential by connecting every cell to its neighbors. To this end, the system is left to evolve for a time , where the characteristic bioelectric time for and . The electrical capacitance C characterizes here the individual cell sensitivity to bioelectrical changes. Then, at the simulation time t = 0 s the potassium channel is blocked, which is modeled by imposing , and the system evolves to different bioelectrical states, first the abnormally depolarized state and then the repolarized state, for long times (Fig. 2).

In Fig. 5, the conductances and the typical parameter values are , k = pol and dep, , , , and , with the threshold potentials for all channels and cells. Here, is the thermal potential, where R is the gas constant, T is the temperature, and F is the Faraday constant16. The pol channel gating charge is zpol = 1 and its conductance (0 for the case of the pol channel blocking). The respective dep channel values are zdep = 2 and . For the res channel, zres = 0.5, with . In the deterministic model, the genetic regulatory networks are characterized by for the mRNA transcription and protein translation rates, respectively, and for the respective degradation rates, with k = pol and dep. For the res channel, (deterministic model) and (stochastic simulations), , and . The res protein reference concentrations are (deterministic model) and (stochastic simulations). The intercellular conductance can take zero, (isolated cells), and non-zero (connected cells) values, with V0 = 10 mV and Vth = 20 mV. Note that we write all conductance values in terms of a reference conductance . The coupling between the bioelectrical and genetic networks of Fig. 5 provides a biophysical scheme based on environmentally-dependent cell states and transitions, including bistability, oscillatory, and synchronization phenomena, that impact on the different gene expression options16.

The evolution of the system is calculated from the local equations of Fig. 5 by using a finite-difference scheme. First, maximum variation for the cell potential at every time step is set at . Similarly, maximum variations of the res channel mRNA and protein concentrations are set at and , respectively. For the sake of the numerical control, we consider the minimum and maximum time steps and , respectively. Then, the time variation of every cell potential Vi is calculated from the equations of Fig. 5 for each simulation time. The time step is (k = pol, dep, res) for all i, keeping . The time, the cell potentials, and the mRNA and protein concentrations are updated as and , and , (k = pol, dep, res), respectively. The calculation proceeds until the desired time is reached. Because typical electrical times can be in the range 0.1–1 s for a reference channel conductance and a cell capacitance C = 10 - 100 pF16, model times of the order of 105 in units correspond to times in the range 3–30 h.

To describe the cell adaptive exploration, we have implemented a stochastic algorithm based on a plausible assumption: the tight control of channel expression should be relaxed under cell stress to permit the adaptation. This algorithm implements and exploratory dynamics based on the stress-driven bioelectrical feedback. The evolution of the systems proceeds similarly as in the deterministic model, i.e. by numerically solving the equations of Fig. 5. However, in this case the transcription rates are changed stochastically every time step after the blocking of the pol channel to simulate the adaptive process. To this end, we introduce the relative potential deviation of cell i as , where is the cell potential and is the target physiological potential, which corresponds to the value before the pol channel blocking. Then, we generate two random numbers per channel following a uniform distribution of values. The first one () is in the range , with , and accounts simply for the basal noise. The second one () is in the range and incorporates a biased noise term because of the small asymmetry provided by the sign (sgn) function. This term results in a biased exploration, according to our assumption that the cell potential should act as a master regulator.

The stochastic algorithm progresses according to the following process. When , the potential of cell i is depolarized with respect to its target value . Then, we change the dep and res channel transcription rates by and , respectively. In this way, the dep channel rate changes following the basal noise while the res channel rate changes following both the basal noise and biased exploration. On the contrary, when , the cell i potential is polarized with respect to its target value. Then, we change the dep and res channel transcription rates by and , respectively. In this case, it is the dep channel rate that changes following both the basal noise and biased exploration while the res channel rate changes following the basal noise only. The small asymmetry in allows the channel rates to change with time in order to compensate for the difference between the cell potential and the target physiological potential. When (depolarized cell) the parameter makes to increase in order to repolarize the cell towards . When (hyperpolarized cell), makes to increase in order to depolarize the cell towards . This scheme also prevents an overshooting of the above repolarization/depolarization processes. To this end, a maximum value of should be introduced to prevent an uncontrolled increase. Note that the changes in the pol channel rate could not modify the cell state, because this channel is blocked.

Wherever the potential of cell i remains deviated from its target potential a value larger than for a time period larger than , the cell is assumed to be non-viable and dies. As could be expected, increasing this viability window results in a smaller number of non-viable cells in the simulations.

We are grateful to the Reviewers for their helpful comments. J. C. and S. M. acknowledge the Ministerio de Ciencia e Innovación (Spain) and the European Regional Development Funds (FEDER), project No PID2022-139953NB-I00. M. L. acknowledges the support by the Templeton World Charity Foundation (TWCF0606) and the Guy Foundation Family Trust 103733-00001.

J.C., M.L. and S.M. made the initial conceptualization. J.C., and S.M. devised the methodology employed. J.C. made the simulations. M.L. and S.M. supervised the research. S.M. wrote the original draft and J.C. prepared the figures. All authors reviewed and edited the manuscript.

The algorithm used here was initially implemented using Python and Javascript as computing languages.For the sake of efficiency and to take advantage of the Template Numerical Toolkit (TNT) and the JAMA/C++ Linear Algebra Package (https://math.nist.gov/tnt/index.html), the code was eventually implemented in C++. The data obtained from the simulations and analysed during the current study are available from the corresponding authors on reasonable request.

M.L. is cofounder of the company Morphoceuticals Inc., which operates in the regenerative medicine and hopes to do regeneration by bioelectric technology. The authors declare no other competing interests.

The algorithm used here was initially implemented using Python and Javascript as computing languages.For the sake of efficiency and to take advantage of the Template Numerical Toolkit (TNT) and the JAMA/C++ Linear Algebra Package (https://math.nist.gov/tnt/index.html), the code was eventually implemented in C++. The data obtained from the simulations and analysed during the current study are available from the corresponding authors on reasonable request.