# Multi-physics modeling for ion homeostasis in multi-compartment plant cells using an energy function

Guillaume Mestdagh, Alexis de Angeli, Christophe Godin

► To cite this version:

Guillaume Mestdagh, Alexis de Angeli, Christophe Godin. Multi-physics modeling for ion homeostasis in multi-compartment plant cells using an energy function. PLoS Computational Biology, 2025, 21 (11), pp.e1013474. 10.1371/journal.pcbi.1013474 . hal-04901993v2

HAL Id: hal-04901993

https://hal.inrae.fr/hal-04901993v2

Submitted on 16 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multi-physics modeling for ion homeostasis in multi-compartment plant cells using an energy function

Guillaume Mestdagh¹, Alexis De Angeli², and Christophe Godin¹

¹Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon1, CNRS, INRAE, Inria, F-69342 Lyon, France

²IPSIM, Univ Montpellier, CNRS, INRAE, Institut Agro, 34060 Montpellier, France

July 16, 2025

## Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>A multi-physics model based on an energy function</td><td>5</td></tr><tr><td>2.1</td><td>Energy function and potentials</td><td>7</td></tr><tr><td>2.2</td><td>Membranes, transporters, directions</td><td>12</td></tr><tr><td>2.3</td><td>Numerical implementation</td><td>16</td></tr><tr><td>3</td><td>Numerical simulations involving a simplified guard cell model</td><td>18</td></tr><tr><td>3.1</td><td>Description of a minimal model representing a guard cell</td><td>19</td></tr><tr><td>3.2</td><td>Numerical simulations involving the minimal guard cell model</td><td>21</td></tr><tr><td>3.3</td><td>Impact of hydrogen buffering mechanisms on simulation results</td><td>26</td></tr><tr><td>4</td><td>Sensitivity analysis using second derivatives of energy</td><td>27</td></tr><tr><td>4.1</td><td>Second derivatives of the energy function: a physical analogy</td><td>27</td></tr><tr><td>4.2</td><td>Comparison between Hessian matrices and simulation results</td><td>30</td></tr><tr><td>4.3</td><td>An interpretation for turgor pressure building in plant cells</td><td>35</td></tr><tr><td>5</td><td>Discussion and Conclusion</td><td>36</td></tr><tr><td>6</td><td>Supporting material</td><td>39</td></tr></table>
Abstract

Plant cells control their volume by regulating the osmotic potential of their cytoplasm and vacuole. Water is attracted into the cell as the result of a cascade of solute exchanges between the cell subcompartments and the cell surroundings, which are governed by chemical, electrostatic and mechanical forces. Due to this multi-physics aspect and to couplings between volume changes and chemical effects, modeling these exchanges remains a challenge that has only been partially addressed. As interest for multi-compartment models grows in the plant cell community, this challenge calls for new modeling strategies. In this paper, we introduce an energy-based approach to couple chemical, electrical and mechanical processes taking place between several subcompartments of a plant cell. The contributions of all physical effects are gathered in an energy function, which allows us to derive the equations satisfied by each variable in a systematic way. We illustrate the properties of this modular, unified approach on the modeling of ion and water transport in a guard cell during stoma opening. We represent the stoma opening process as a quasi-static evolution driven by hydrogen pumps in the plasma and vacuolar membranes, and we show that the new formalism explains why the system varies in a particular direction in response to perturbations. Additional numerical simulations allow us to investigate the role of each hydrogen pump in this process. Altogether, we show that this energy-based approach highlights a hierarchy between the forces involved in the system, and to dissect the role of each physical effect in the complex behavior of the system.

**Keywords** Cell turgidity regulation, ion transport, proton pump, cell physics, energy-based formalism.

Author summary

Osmosis is a physical effect by which water is attracted into compartments with a high concentration of solute. Plant cells exploit osmosis to attract water into their subcellular compartments and increase their sizes. They use complex mechanisms, involving the circulation of several ion species between compartments, governed by chemical and physical phenomena of multiple natures. Due to the complexity of their interactions, these mechanisms are difficult to study experimentally, and require the development of quantitative modeling approaches. In this paper, we propose a theoretical model whose goal is to show how these ions exchanges result in water finally flowing into the cell. Particular attention is paid to the competition between the various physical phenomena, their priorities and the way each of them influences the system.
# 1 Introduction

To achieve key living functions, plant cells change their size and shape by regulating the amount of water entering their vacuole [1, 2, 3, 4]. For this, they actively control the osmotic potential of their cytoplasm and vacuole by importing or exporting ions or sugars across the membranes delimiting these compartments (the plasma and vacuolar membranes, respectively). Solutes are transported by a wide variety of specialized proteins located on each membrane, including active pumps, channels, exchangers, etc. Due to their specific stoichiometry, these transport mechanisms define constraints on possible transfers of species between subcellular compartments. The forces that drive solute transfer are of multiple natures. They integrate chemical, electrical, and mechanical components, whose combination determines the dynamics of this multi-physics, multi-membrane system.

Modeling approaches have been widely used to help analyze the behavior of such complex systems, as a complement to experimentation. A first category of approaches focuses on electro-chemical exchanges across membranes using ordinary differential equations. In these models, the variation of ion concentrations in each cell compartment is defined by the ion fluxes crossing membranes through a range of transport proteins. Vitali et al. [5], for example, use simulation to identify the most significant parameters controlling the volume dynamics of beetroot vacuoles. In a recent work, Li et al. [6] account for the multimembrane structure of the cell, and simulate the ion fluxes across both the plasma and vacuole membranes of a guard cell in response to variations in external potassium and chloride concentrations. Note that similar studies have also been carried out to model ion transport in animal or fungi cells [7, 8]. Many cell physiology models use a formalism inspired by chemical kinetics to represent the dynamics of each transport process. Gerber et al. [9] opt for a different formalism inspired by thermodynamics, introduced in the 1950s [10, 11, 12], where the flux across a membrane is proportional to the electrochemical potential gradient between the two connected compartments. They use a parameter estimation procedure to fit their simulation results to experimental measures. The variety of modeling approaches used to determine the flux of ions between compartments illustrates the fact that the modeling of transport dynamics in cells remains a challenging and discussed topic.

A second category of approaches aims at modeling the mechanical response of the cell wall to an increasing turgor pressure inside the cell. In the case of cell growth, this response is irreversible, as volume increase triggers wall remodeling and cell division. Complex mechanical models proposed in this context are reviewed in [13, 14]. These models account for the influence of the composite structure of cell walls [15], membrane integrity [16] or enzyme activity [17]. Mechanical response can also be reversible, as in the case of guard cells. Located on the surface of plant leaves and stems, guard cells work in pairs to control the aperture of small pores called stomata. These specialized cells are able to repeatedly get inflated and deflated to open and close stomata, and regulate plant transpiration and carbon dioxide absorption [18, 19]. Many mechanical studies examine the elastic properties of walls in guard cells [see the review by 20]. Meckel et al. [21] discuss whether guard cells widen or elongate as their volume increases. By comparing finite element simulations and experimental measures, Woolfenden et al. [22] estimate the stiffness parameter distribution of the cell wall. Similar approaches are followed in [23, 24]. Carter et al. [25] study the role of wall stiffening at the extremities of guard cells in the
stoma opening process, while Jaafar et al. [26] focus on the role of cell wall anisotropy.

The models mentioned above either describe electrochemical or mechanical phenomena. On the one hand, ion exchange models involve an osmotic pressure but rarely describe water movements between compartments, whereas on the other hand, mechanical models predict cell deformation as a result of a turgor pressure whose origin is not specified. The interplay between these two categories of effects is often not taken into account, as it represents an additional level of complexity. Also, the multi-compartment aspect of plant cells is rarely considered. Recent studies [27, 6] highlight the need to include it in modeling approaches, to better understand how cells synchronize the regulation of their membranes. This new trend in plant cell physiology calls for models and formalisms able to provide novel insight into multi-compartment systems.

The challenge of building multi-membrane models that couple electrochemistry with mechanics is addressed in the OnGuard model [28]. This well-established framework aims to provide realistic simulations of solute exchanges in a guard cell. To connect the guard cell physiology with the stoma aperture, the authors use a linear pressure-volume law, along with a linear volume-aperture relation, whose coefficients are empirically set based on experimental results. OnGuard has been implemented as a software platform, that has proven useful to reproduce experimental results [29], and to drive new experiments by predicting guard cell behavior [30, 31, 32]. To run a simulation, the user chooses among an exhaustive list of membrane transport proteins and sets the dynamic parameters of each reaction law. Results are available through different plots showing the time evolution of the guard cell volume and aperture, osmolyte concentrations and fluxes, among many other physical quantities.

In this paper, similarly to OnGuard, our main objective is to develop a multi-membrane, multi-physics model for the regulation of plant cell volume control. However, we follow a different and complementary approach. While OnGuard focuses on exhaustiveness and resemblance to experimental results, we center our attention on simplicity and modularity. We seek to integrate electrochemical and hydro-mechanical processes in an explicit and mechanistic way, while relying on few parameters, whose relative effect on the simulations can be qualitatively and quantitatively assessed. To meet these ambitions, we unify the various physical processes by gathering their contributions into a common global energy function. This energy function naturally brings some modularity to the model, as its components can be added or removed without impact on the other components of the system. In addition, it keeps generic features (e.g. implementation of the main physical processes) well separated from specific ones (e.g. choice of specific transport stoichiometry, etc.). We exploit this approach to construct a minimal model representing a guard cell during stoma opening, and simulate ion transfer between its subcellular compartments. To keep the model simple, we assume that the underlying physical processes follow at every moment a moving equilibrium, which progressively evolves in response to successive perturbations. This makes it possible to neglect the dynamics of ion transport and to keep fewer parameters in the model. In such quasi-static simulations, time-dependent ordinary differential equations are replaced with a series of nonlinear equations that are solved to compute the sequence of equilibrium states. Altogether, the use of an energy landscape and quasi-static simulation allows to obtain a simple geometrical interpretation of the system evolution. We will show that the resulting formulation is flexible and easily provides physical insight into the role of the various physical processes involved in the system regulation.
## 2 A multi-physics model based on an energy function

We now introduce the formalism used to construct a simple representation of a plant cell. We reduce the cell to a set of nested compartments connected by transporters. For clarity reasons, we use the term transporter to designate any process that transports solute or water across a membrane. This definition includes pumps, channel, coupled ion transport, and even water permeating across membranes. As mentioned in Section 1, we seek to compute a sequence of equilibrium states for the system, that is parametrized by the chemical extent of some controlled transport reactions. In other words, a resulting equilibrium state for the whole system is associated with a given state of the active transport processes. As active transporters work, the system equilibrium state constantly adjusts in consequence.

From now on, we reduce an object as complex as a plant cell to a multi-membrane complex, which is an arrangement of one, two, or possibly more nested compartments enclosed by membranes. Each compartment contains a certain chemical amount (in mol) of chemical species such as ions and water. In this article, water can circulate across membranes, while the transport of ions between compartments is allowed by several active and passive transporters. The multi-membrane complex is surrounded by an external environment with fixed ion concentrations, which acts as a reservoir. A plant cell is typically represented by a complex with two compartments (the cytoplasm and the vacuole). For instance, the system shown in Figure 1, involves two membranes, four chemical species (hydrogen, chloride, potassium, water) and eight transporters.

The state of the system at a given moment is entirely defined by the chemical amount of each considered chemical species in each compartment. Of course, to perform simulations involving a multi-membrane model, some fixed parameters need to be defined, such as the initial content of all compartments or initial electric potentials across membranes, but these parameters are not variables, and they remain fixed along the simulation. The chemical amount of species A in compartment i is denoted by $n_{i,A}$, and all chemical amounts are gathered in the state vector $\mathbf{n} \in \mathbb{R}^N$. Here, the size N of the state vector is the product of the number of chemical species with the number of compartments. For instance, the state vector for the model from Figure 1, of size N = 8, reads

$$ \mathbf{n} = (n_{1,H^+} \ n_{1,Cl^-} \ n_{1,K^+} \ n_{1,H_2O} \ n_{2,H^+} \ n_{2,Cl^-} \ n_{2,K^+} \ n_{2,H_2O}) \in \mathbb{R}^8. \quad (2.1) $$

Note that the variables in the state vector are first ordered by compartment, and then by chemical species, and that the components correspond to amounts of chemical species rather than to concentrations.

An equilibrium state is defined as the state that minimizes an energy function among a range of states attainable by the system. More formally, the equilibrium state is the solution to a constrained optimization problem of the form

$$ \text{minimize } G(\mathbf{n}) \quad \text{subject to constraint} \quad \mathbf{n} \in \mathcal{V}_{\mathbf{x}}, \quad (2.2) $$

where G is an energy function to minimize, n is the variable that describes the system state, and $\mathcal{V}_{\mathbf{x}}$ is the range of states attainable by the system, which depends on an input parameter x.
<figure><img src="image_2.png" /><figcaption>Figure 1: Representation of a guard cell as a multi-membrane complex. The model involves two membranes and a hand-picked selection of transporters. Only the outer membrane offers mechanical resistance to volume changes, which results in an equal pressure in both compartments. Each compartment contains a certain amount of hydrogen, chloride, potassium and water. All these chemical amounts compose the system state, from which the energy function is evaluated. (Gear image adapted from Issi at openclipart.org).</figcaption></figure>
Below, we give details about the various components of this minimization problem, starting with the energy function.

## 2.1 Energy function and potentials

To understand the dynamics of this physical system, we seek to describe how its protagonists store, exchange or dissipate energy. In particular, the energy that is stored and does not correspond to a current motion is called potential energy. It is the energy paid by the system just to be in a certain state. The potential energy function of a system depends on the system state only, and it describes the forces that influence the system. These forces make the physical system spontaneously evolve to decrease its potential energy, generally transforming it into heat. A state where the potential energy is at a local minimum is called an equilibrium state. When a physical system has reached an equilibrium, it ceases to evolve, and when it is taken away from this equilibrium, it tends to go back to this equilibrium (we will only meet stable equilibria in this study). From now on, we simply use *energy function* to refer to the potential energy function. The potential energy function should not be confused with potentials (chemical potential, electric potential), that rather correspond to derivatives of the energy function. A chemical potential is expressed in J/mol, an electric potential in J/C and a pressure in J/m$^3$, but the potential energy is expressed in J.

Defining an energy function for a physical system is a very convenient modeling practice, as energy is a common currency between a wide range of physical processes. In general, the system energy function is the sum of several terms, with each term corresponding to a specific physical process. In this paper, the energy function, is split into three terms, as we consider three types of physical effects. First, chemical effects tend to equilibrate the concentration of solutes between both side of each membrane. Second, electrostatic effects oppose to the building of large electric potential variations across membranes. Finally, mechanical effects act on volume changes.

Though the energetic formalism is compatible with several choices of shape for the cell, geometry plays a significant role when it comes to deriving an actual expression for energy function terms. In this paper, we stick to a very simple geometry for each membrane, namely a cylinder that elongates in one direction only, along its revolution axis, so that one-dimensional considerations are sufficient to define the elastic deformation energy. When the cell volume changes, only the cylinder lateral surface is subject to strain in the direction of the elongation, and the cylinder volume and lateral surface area remain proportional to its length. We denote by $V_i(\mathbf{n})$ the volume enclosed in the $i$-th membrane, and by $A_i(\mathbf{n})$ the area of its lateral surface. For a cylinder with fixed radius $r$, the lateral surface area and the volume satisfy $A_i(\mathbf{n}) = 2V_i(\mathbf{n})/r$. Note that the $i$-th membrane encloses compartment $i$, but also all compartments with index $j > i$.

We now give a little more details about each term of the energy function.

**Chemical energy** The chemical energy term measures a discrepancy between the system state and a chemical equilibrium. Chemical forces usually determine in which direction a chemical reaction occurs [33, Chapter 8]. Whereas we do not consider any chemical transformation in this model, chemical forces still play a role in determining which way ions circulate through
a transporter. Namely, the chemical energy takes larger values when the concentration of a species is very different between two compartments, and chemical forces tend to promote, for each species, a uniform concentration in the compartments. Chemical forces are responsible for the osmosis process, i.e. attracting water in compartments where the solute concentration is high.

We use a standard expression for the chemical energy term, based on each chemical potential $\mu_{i,A}(\mathbf{n})$ associated with species A in compartment i. When chemical transformations occur, the chemical potential determines how the chemical energy changes when species A is produced or consumed in compartment i, and its classical expression reads

$$ \hat{\mu}_{i,A}(\mathbf{n}) = \mu_{A}^{\circ} + RT \ln(\gamma_{A}y_{i,A}(\mathbf{n})), \quad \text{with} \quad y_{i,A}(\mathbf{n}) = \frac{n_{i,A}}{\sum_{R} n_{i,R}}, \qquad (2.3) $$

where $y_{i,A}(\mathbf{n})$ the mole fraction of species A in compartment i, and $\mu_{A}^{\circ}$ and $\gamma_{A}$ represent the chemical potential of pure species and the activity coefficient of species A, respectively [34, Section 6.3]. However, we only consider solute transport in this study, and an increase of $n_{i,A}$ without variation of other quantities means that particles of A have been taken from the external environment and transported into compartment i. The corresponding marginal variation rate of chemical energy involves the chemical potential of A in the external environment. It reads

$$ \mu_{i,A}(\mathbf{n}) = \hat{\mu}_{i,A}(\mathbf{n}) - \hat{\mu}_{\text{ext},A} = RT \ln \left( \frac{y_{i,A}(\mathbf{n})}{y_{\text{ext},A}} \right). \qquad (2.4) $$

Here, the external environment is considered as a reservoir, where each species A is present with the constant mole fraction $y_{\text{ext},A}$, so that the chemical potential $\hat{\mu}_{\text{ext},A}$ of A in the external environment is constant. The expression for chemical energy is chosen so that $\frac{\partial G_{\text{chem}}}{\partial n_{i,A}}(\mathbf{n}) = \mu_{i,A}(\mathbf{n})$, which yields

$$ G_{\text{chem}}(\mathbf{n}) = \sum_{i} \left( \sum_{A} n_{i,A} RT \ln \left( \frac{y_{i,A}(\mathbf{n})}{y_{\text{ext},A}} \right) \right) = \sum_{i} \left( \sum_{A} n_{i,A} \mu_{i,A}(\mathbf{n}) \right). \qquad (2.5) $$

Note that $\mu_{i,A}$ only depends on chemical amounts in compartment i, so that (2.5) is actually made of independent terms, with each term depending on mole fractions in one compartment. From expressions (2.5) and (2.4), we can check that the chemical energy function is at a minimum when all chemical potentials are zero, i.e. mole fractions in each compartment are equal to external mole fractions.

**Electrostatic energy** The electrostatic term of the energy function aims to penalize large differences of electric potential across each membrane. When a charged particle crosses a membrane from one compartment to another, the net transfer of electric charge results in a charge imbalance across the membrane and a difference of electric potential between the two compartments. This phenomenon, also known as membrane polarization, tends to oppose to the passage of more similar charges in the same direction. From a modeling point of view, a membrane can be considered as an electric capacitor, i.e. a thin insulating layer between two charged compartments, allowing charges to accumulate along the membrane when a difference of potentials
<figure><img src="image_3.png" /><figcaption>Figure 2: Electric representation of a two-membrane complex. Each capacitor represents a membrane and each compartment corresponds to a node in the circuit. Variations of ion amounts in each compartment result in adding or removing charges in these compartments. The circuit representation suggests that adding charged particles into compartment 2 contributes to charging both capacitors $C_1$ and $C_2$. This circuit analogy is used to evaluate the electrostatic energy term.</figcaption></figure>

occurs. As a consequence, the electrostatic energy value is the sum of the energy of all capacitors. The electrical capacitance of a membrane is proportional to its area. The capacitance per unit membrane area in biological cells is typically 1 µF/cm² [35, Section 2.6].

The electric structure of a multi-membrane complex can be represented by an equivalent circuit, shown in Figure 2 in the two-membrane case. Each compartment corresponds to a node in the circuit, while the external environment is represented by the ground. Compartments are connected by capacitors that represent the membranes, while the two current sources symbolize the action of adding charges from the external environment into one of the compartments. To understand the relationship between the content of each compartment and the state of capacitors, let us apply Kirchhoff’s circuit laws: if we add a charge $dq_1$ into compartment 1 and $dq_2$ into compartment 2, then the charge stored in capacitor $C_2$ changes by $dq_2$, while the charge stored in capacitor $C_1$ changes by $dq_1 + dq_2$ (here, charge displacements act like current intensities). In other words, the charge $Q_i(\mathbf{n})$ stored in the capacitor $C_i$ depends on the cumulated electric charge in all compartments enclosed inside the i-th membrane. In particular, transporting charged particles from compartment 1 to compartment 2 corresponds to the situation $dq_1 = -dq_2$, which does not change the charge stored in capacitor $C_1$.

Let us now derive an expression for the electrostatic energy stored in capacitor $C_i$. If we denote by $\mathbf{n}_k = (n_{k,A}, n_{k,B}, \dots)$ the vector of all chemical amounts in compartment $k$, the charge contained in compartment $k$ is defined by

$$q_k(\mathbf{n}) = q_k^0 + F\mathbf{z}^T(\mathbf{n}_k - \mathbf{n}_k^0), \qquad (2.6)$$

where $\mathbf{z} = (z_A, z_B, \dots)$ stores the charge of each species, and $F$ is the Faraday constant. The
parameters $q_k^0$ and $\mathbf{n}^0$ are fixed as part of the initial conditions. The expression for the charge stored in capacitor $C_i$ reads

$$Q_i(\mathbf{n}) = \sum_{k \ge i} q_k(\mathbf{n}) = Q_i^0 + \sum_{k \ge i} F \mathbf{z}^T (\mathbf{n}_k - \mathbf{n}_k^0), \qquad (2.7)$$

where $Q_i^0$ denotes the stored charge in initial conditions. The energy stored in capacitor $C_i$ is defined using the classic formula

$$G_{i,elec}(\mathbf{n}) = \frac{1}{2} \frac{Q_i(\mathbf{n})^2}{C_i(\mathbf{n})} = \frac{1}{2C_i(\mathbf{n})} \left( Q_i^0 + \sum_{k \ge i} F \mathbf{z}^T (\mathbf{n}_k - \mathbf{n}_k^0) \right)^2. \qquad (2.8)$$

Note that the initial membrane charge $Q_i^0$ can be deduced from experimental measures of initial membrane potentials. The expression for the capacitance $C_i(\mathbf{n})$ depends on the cell geometry. As mentioned above, the cell is modeled as a cylinder, and for simplicity we only consider its lateral surface to evaluate the membrane capacitance, which is defined by

$$C_i(\mathbf{n}) = C_m A_i(\mathbf{n}) = 2 \frac{C_m}{r} V_i(\mathbf{n}), \qquad (2.9)$$

where $C_m$ is the capacitance per unit surface.

Finally, we evaluate the total electrostatic energy by adding up the energy stored in each capacitor, i.e.

$$G_{elec}(\mathbf{n}) = \sum_{i} G_{i,elec}(\mathbf{n}). \qquad (2.10)$$

**Mechanical energy and cell geometry** The mechanical energy term represents the elastic energy stored in the cell wall as the cell volume increases. While cell growth is an irreversible process involving cell wall elasticity and plasticity at the same time, guard cell deformation during stoma opening is reversible, as the stoma opening/closing cycle happens several times, at least once a day. As we focus in this paper on the quasi-static modeling of a guard cell, the only mechanical effect we consider is elasticity.

The mechanical energy stored in a membrane is a function of the volume it encloses. When subject to a volume change, the membrane adopts a resulting mechanical configuration, which corresponds to a value for its mechanical energy. In this paper, the only contributor to compartment volume is water, so that the volume enclosed in the i-th membrane is defined by

$$V_i(\mathbf{n}) = v_{H_2O} \sum_{k \ge i} n_{k,H_2O}, \qquad (2.11)$$

where $v_{H_2O}$ is the molar volume of water. Keep in mind that $V_i(\mathbf{n})$ is not the volume of compartment i, but the cumulated volume of all compartments enclosed in the i-th membrane. Also, in the context of plant cells, we assume that the only membrane that offers mechanical resistance to volume changes is the outer membrane (membrane 1), as in plant cells the plasma membrane is coupled with the cell wall, while the vacuole membrane (membrane 2) is only a lipid bilayer. The resulting elastic deformation energy only depends on the whole cell volume $V_1(\mathbf{n})$, which is
the sum of all compartment volumes. In other words, defining the mechanical energy term boils down to defining a function of one variable, sometimes referred to as a pressure/volume law. In terms of pressure, it means that no membrane except the first one applies mechanical forces on the fluid it contains, and the mechanical pressure is the same in every compartment. If we define the pressure $P = \partial G_{\text{mech}} / \partial V_1$, by looking at (2.11) we obtain

$$ \frac{\partial G_{\text{mech}}}{\partial n_{1,H_2O}}(\mathbf{n}) = \dots = \frac{\partial G_{\text{mech}}}{\partial n_{N_{\text{comp}},H_2O}}(\mathbf{n}) = v_{H_2O} P(\mathbf{n}), \qquad (2.12) $$

which means that it has the same cost in mechanical energy to add some water into any compartment.

Without choosing a cell geometry, it is difficult to be more specific about the mechanical energy term. We now exploit the cylindrical cell shape to come up with an expression. When the cell expands, the strain in the lateral surface is given by the expression $\epsilon(\mathbf{n}) = V_1(\mathbf{n}) / V_1(\mathbf{n}^0) - 1$. We denote by $\Sigma^0$ the initial cylinder lateral surface and by $A_i(\mathbf{n}^0)$ its area. Using a linear constitutive law, the energy integrated over the whole lateral surface reads

$$ G_{\text{mech}}(\mathbf{n}) = \int_{\Sigma^0} \frac{1}{2} E_s \epsilon(\mathbf{n})^2 \, ds = A_1(\mathbf{n}^0) \times \frac{E_s}{2} \left( \frac{V_1(\mathbf{n})}{V_1(\mathbf{n}^0)} - 1 \right)^2, \qquad (2.13) $$

where $E_s$ is a surface elastic modulus (in Pa·m).

Though we used a simple one-dimensional model for the cell geometry, the energetic formalism is completely compatible with alternative descriptions of the cell expansion [21] or more complex geometries [22, 20].

**Enforcing equilibrium in starting conditions** To ensure consistency between the model and experimental conditions, we need to start a simulation from any user-specified initial state $\mathbf{n}_0$, which can be a state measured at the beginning of an experiment. However, in a quasi-static model, $\mathbf{n}_0$ must be an equilibrium state of the system, i.e. it should correspond to a minimum of the energy function $G$. To make sure that $\mathbf{n}_0$ represents an equilibrium state, we add a linear correction term $-\mathbf{f}_0^T \mathbf{n}$ to the energy function $G$, where $\mathbf{f}_0$ is adjusted so that $\nabla G(\mathbf{n}_0) = 0$, which means that $\mathbf{n}_0$ is a minimum of $G$. Namely,

$$ \mathbf{f}_0 = \nabla G_{\text{chem}}(\mathbf{n}_0) + \nabla G_{\text{elec}}(\mathbf{n}_0) + \nabla G_{\text{mech}}(\mathbf{n}_0). \qquad (2.14) $$

The force $\mathbf{f}_0$ can be interpreted as the contribution of the environment to the energy function. A plant cell is the theater of many processes, associated with metabolism or other functions that are not explicitly modeled. From the point of view of our model, we assume that all this activity can be summarized by the constant force $\mathbf{f}_0$. For a given $\mathbf{n}_0$, the total energy function that is actually minimized during the simulation is defined by

$$ G(\mathbf{n}) = G_{\text{chem}}(\mathbf{n}) + G_{\text{elec}}(\mathbf{n}) + G_{\text{mech}}(\mathbf{n}) - \mathbf{f}_0^T \mathbf{n}. \qquad (2.15) $$

Though this energy function depends on $\mathbf{n}_0$, we keep calling it $G$ in the remaining of this section, as $\mathbf{n}_0$ is a fixed parameter of the model.
<figure><img src="image_4.png" /><figcaption>Figure 3: A toy multi-membrane complex with one membrane and three transporters. The hydrogen pump is an active transporter, while the hydrogen/chloride symporter is passive. Water transport across the membrane is also considered passive. As there are three species and one compartment, the system state is a vector of size 3. We use this toy problem to illustrate the proposed quasi-static setting.</figcaption></figure>

## 2.2 Membranes, transporters, directions

The energy function described in Section 2.1 defines the forces that influence the system. It is completely agnostic of the disposition of transporters. If the system was free to reach any state $n \in \mathbb{R}^N$, it would converge toward the global minimum of the energy function $G$. However, some states cannot be reached by the system, as chemicals can only circulate between compartments through transporters, which exhibit a specific stoichiometry. The disposition of transporters defines the set of states attainable by the system, denoted by $\mathcal{Y}_x$ in (2.2).

In this subsection, we describe the connection between the disposition of transporters and the set of attainable solutions. For the sake of simplicity, we use a simple one-membrane complex to illustrate the chosen formalism. This toy system, shown in Figure 3, involves one compartment, three chemical species (hydrogen, chloride, water) and three transporters. The same formalism will later be extended, *mutatis mutandis*, to a complex with several membranes.

**Points and directions in the state space** In the one-membrane complex, the system state is defined by the chemical amount of each species A in the single compartment. We use the slight abuse of notation $n_A = n_{1,A}$, so that the state vector now reads

$$ \mathbf{n} = (n_{H^+} \quad n_{Cl^-} \quad n_{H_2O}) \in \mathbb{R}^3. \qquad (2.16) $$
The vector **n** may be interpreted as a point in a three-dimensional state space, whose coordinates are the amounts of considered chemical species. Also, a variation of chemical amounts corresponds to a displacement along a certain direction in this space, represented by a vector of size 3. For instance, in Figure 3, the first transporter is a hydrogen pump. It ejects hydrogen from the compartment to the external environment, which results in the hydrogen amount $n_{H^+}$ decreasing. As $n_{H^+}$ is the first coordinate of the state vector **n**, this is represented in the state space by **n** evolving along the vector

$$ \mathbf{a} = (-1 \ 0 \ 0), \qquad (2.17) $$

where we have used the same ordering as (2.16) for chemical species. Once the hydrogen pump has ejected $x$ moles of hydrogen (and no other transport occurred), the new state now reads

$$ \mathbf{n} + \mathbf{a}\mathbf{x} = (n_{H^+} - x \ n_{Cl^-} \ n_{H_2O}). \qquad (2.18) $$

The second transporter is called a symporter because it transports two species in the same direction, at the rate of one hydrogen ion for two chloride ions, while the third transporter symbolizes water permeating across the membrane. In a similar fashion as previously, each one of these two transporters allows the state to evolve in a particular direction. These directions, denoted by $\mathbf{v}_1$ and $\mathbf{v}_2$, are defined after the stoichiometry of the two transporters (see Figure 3), namely

$$ 1 \mathbf{H}^+ + 2 \mathbf{Cl}^- + 0 \mathbf{H}_2\mathbf{O} \leftrightarrow \mathbf{v}_1 = (1 \ 2 \ 0) \qquad (2.19) $$

$$ \text{and} \quad 0 \mathbf{H}^+ + 0 \mathbf{Cl}^- + 1 \mathbf{H}_2\mathbf{O} \leftrightarrow \mathbf{v}_2 = (0 \ 0 \ 1). \qquad (2.20) $$

In a general multi-membrane complex, the state can only evolve in directions defined by the stoichiometric coefficients of transporters, which restrict the range of attainable states. For instance, in the toy model, if we block the hydrogen pump, then the amount of chloride cannot change independently of hydrogen, as the only way out for chloride is through the hydrogen/chloride transporter. In particular, the value $2n_{H^+} - n_{Cl^-}$ must remain constant. Here, the state **n** can only move along two directions $\mathbf{v}_1$ and $\mathbf{v}_2$, and therefore the set of attainable states is a subspace of dimension two.

**Constrained energy minimization** The hydrogen pump is a particular transporter: whereas the other transporters just passively allow ions and water to cross the membrane at the whim of electrochemical potentials, the pump consumes some energy from an external source to eject hydrogen from the compartment, sometimes against chemical potentials. We call the hydrogen pump an *active* transporter while other transporters are *passive*. This distinction is critical for the remaining of this section, as the directions associated with passive and active transporters play different roles in the definition of $\mathcal{Y}_x$. The stoichiometry of passive transporters defines the directions that span $\mathcal{Y}_x$, while active transporters define the displacement of $\mathcal{Y}_x$ in directions orthogonal to itself.

Let us assume for a moment that the active pump is not working, which means that hydrogen and chloride can only circulate together through the hydrogen/chloride transporter. Starting from an initial state $\mathbf{n}^0$, the system evolves toward states of lower energy, but, as mentioned above, the state vector **n** can only move along the two linearly independent directions $\mathbf{v}_1$ and
$v_2$. In other words, $\mathbf{n}$ is restricted to an affine plane $\mathcal{V}_0 \subset \mathbb{R}^3$ containing $\mathbf{n}^0$. If we denote by $\mathbf{w}_1 = \mathbf{v}_1 \times \mathbf{v}_2 = (2, -1, 0)$ a normal vector to $\mathcal{V}_0$, the subspace of attainable states is characterized by the equation

$$ \mathbf{n} \in \mathcal{V}_0 \Leftrightarrow \mathbf{w}_1^T (\mathbf{n} - \mathbf{n}^0) = 0 \Leftrightarrow 2n_{H^+} - n_{Cl^-} = 2n_{H^+}^0 - n_{Cl^-}^0, \quad (2.21) $$

where $\mathbf{w}_1^T$ is the transpose of $\mathbf{w}_1$. As a consequence, when the active transporter does not work, the equilibrium attained by the system is the minimizer of the energy function $G$ among the subspace $\mathcal{V}_0$. Now, assume that the hydrogen pump ejects $x$ moles of hydrogen before stopping, bringing the system to the state $\mathbf{n}^0 + \mathbf{a}x$ (see 2.18). To find an equilibrium, the state $\mathbf{n}$ is now restricted to another subspace, denoted by $\mathcal{V}_x$, whose equation reads

$$ \mathbf{n} \in \mathcal{V}_x \Leftrightarrow \mathbf{w}_1^T (\mathbf{n} - \mathbf{n}^0 - \mathbf{a}x) = 0 \Leftrightarrow 2n_{H^+} - n_{Cl^-} = 2(n_{H^+}^0 - x) - n_{Cl^-}^0. \quad (2.22) $$

It is notable that, for $x \neq x'$, the subspaces $\mathcal{V}_x$ and $\mathcal{V}'_x$ are disjoint ($2n_{H^+} - n_{Cl^-}$ takes different values), and as a consequence, the resulting system equilibria are distinct.

Figure 4 shows a slice view of the space of states in the $(n_{H^+}, n_{Cl^-})$ plane, with representations of the subspaces $\mathcal{V}_x$ corresponding to several values of $x$. The lines in the background represent the level curves of the energy function $G$. The point that minimizes $G$ in a given subspace $\mathcal{V}_x$ is the point where a level curve is tangent to $\mathcal{V}_x$. As $\mathcal{V}_x$ moves, the equilibrium state $\mathbf{n}^*(x)$ is pushed toward higher levels of energy, following the red curved trajectory.

In the toy model from Figure 3, the extent of active transport $x$ is the input variable, while the output is the corresponding equilibrium state $\mathbf{n}^*(x)$. When the amount of hydrogen ejected by the pump describes an interval $[0, x]$, it progressively displaces the passively attainable subspace $\mathcal{V}_x$, and the system equilibrium is continuously brought from $\mathbf{n}^*(0) \in \mathcal{V}_0$ to $\mathbf{n}^*(x) \in \mathcal{V}_x$. A more general multi-membrane complex might feature several active transporters, whose chemical extents are stored in a vector $\mathbf{x} = (x_1, x_2, \dots)$, and whose corresponding directions are stored as the columns of a matrix $\mathbf{A} = (\mathbf{a}_1|\mathbf{a}_2|\dots)$. In that case, the system equilibrium is a function of the current extent of all active transporters, defined by

$$ \mathbf{n}^*(\mathbf{x}) = \arg\min_{\mathbf{n}} G(\mathbf{n}) \quad \text{subject to constraint} \quad \mathbf{n} \in \mathcal{V}_x. \quad (2.23) $$

The space $\mathcal{V}_x$ can always be characterized by an equation of the form

$$ \mathbf{n} \in \mathcal{V}_x \Leftrightarrow \mathbf{W}^T \mathbf{n} = \mathbf{W}^T (\mathbf{n}^0 + \mathbf{A} \mathbf{x}), \quad (2.24) $$

where the matrix $\mathbf{W}$ has $N$ rows and $N - \dim(\mathcal{V}_x)$ columns, and the columns of $\mathbf{W}$ form a basis of $\mathcal{V}_x^\perp$.

**A quasi-static simulation scenario** In this paper, the objective of the simulation is to evaluate the function $\mathbf{x} \mapsto \mathbf{n}^*(\mathbf{x})$ when $\mathbf{x}$ describes a known range of input values. The expected result is an evolution of the quasi-static equilibrium state as active transporters work. In practice, when the input variable $\mathbf{x}$ is multidimensional, we make it follow a one-dimensional trajectory, so that the quasi-static equilibrium follows a one-dimensional trajectory in the space of states. Thus,
<figure><img src="image_5.png" /><figcaption>Figure 4: Slice view of the space of states in the plane ($n_{H^+}$, $n_{Cl^-}$) for the toy problem. The level curves in the background represent the energy function $G$. Once the hydrogen pump has brought the system from $\mathbf{n}_0$ to $\mathbf{n}_0 + x\mathbf{a}$, the passive evolution toward an equilibrium is only possible along the affine subspace $\mathcal{V}_x$, which is directed by $\mathbf{v}_1$ and $\mathbf{v}_2$ ($\mathbf{v}_2$ is out of the plane). The minimum of $G$ in $\mathcal{V}_x$ corresponds to points where the level curves of $G$ are tangent to $\mathcal{V}_x$. As $x$ increases, the plane $\mathcal{V}_x$ continuously moves, dragging the equilibrium state $\mathbf{n}^*(x)$, which follows the thick curved trajectory. Thus, the system state always minimizes the same function $G$, but in a subspace that moves as $x$ increases.</figcaption></figure>
simulation results look like an evolution as a function of a scalar variable. For instance, the example from Figure 1 involves two active transporters. If we run the simulation for $x(s) = (s, s/2)$ with $s$ describing the interval $[0, s_{max}]$, it means that the second pump works twice as slow as the first one. However, keep in mind that $s$ does not represent time here, as our model does not include dynamics. By changing the way $x$ depends on $s$, one can explore several possibilities of synchronization between active transporters.

Choosing a quasi-static scenario is a strong choice, as the individual dynamics of transporters is not taken into account, and transient regimes are not modeled. It corresponds to a limit case where the dynamics of passive transporters and mechanical deformations is very fast compared with the time span of the considered process, so that the system immediately goes back to an equilibrium when perturbed by active transport. Actually, assuming a quasi-static evolution boils down to separating transport processes into two categories: either the transport is slow, and it is not included in the model, or it is sufficiently fast, and it is considered instantaneous. However, our goal here is not to depict a time-dependent evolution. As the dynamic properties of transporters is not always well known, a quasi-static approach is a good way to highlight phenomena that do not depend on the transporter dynamics. Elastic deformations of the cell wall are also considered instantaneous, similarly to the classic Lockhart-Ortega model [36]. Finally, the main advantage of choosing a quasi-static model is that it requires few parameters. While a dynamic simulation requires a function to determine the dynamics of each transporter, sometimes requiring to set many dynamics-related coefficients, the quasi-static model only requires the transporters’ stoichiometry, and the distinction between active and passive transporters.

## 2.3 Numerical implementation

Our simulation code<fnref n="1" /> is written in Python, and we leverage the automatic differentiation possibilities provided by the Jax library [37] to access the derivatives of $G$ at least cost in terms of development. Using automatic differentiation is very convenient when working with variational problems, as resolution algorithms often require to evaluate the higher-order derivatives of energy or cost functions. Here, we were able to run simulations and perform sensitivity analysis with only a limited amount of calculations.

First, let us derive the nonlinear system that is actually solved in the numerical procedure. For a fixed input variable $x$, the corresponding equilibrium is defined by

$$ \mathbf{n}^*(\mathbf{x}) = \arg\min_{\mathbf{n}} G(\mathbf{n}) \quad \text{subject to constraint} \quad \mathbf{W}^\mathrm{T}\mathbf{n} = \mathbf{W}^\mathrm{T}(\mathbf{n}^0 + \mathbf{A}\mathbf{x}), \qquad (2.25) $$

with the notation from (2.23) and (2.24). It can be shown that the solution to (2.25) satisfies some necessary first-order optimality conditions, also known as Karush-Kuhn-Tucker conditions [38, Section 12.3]. These conditions read

$$ \left\{ \begin{array}{l} \nabla G(\mathbf{n}) + \mathbf{W}\lambda = 0 \\ \mathbf{W}^\mathrm{T}(\mathbf{n} - \mathbf{n}^0) = \mathbf{W}^\mathrm{T}\mathbf{A}\mathbf{x}, \end{array} \right. \qquad (2.26) $$

<footnote n="1">¹Code available on github.com/gmestdagh/quasistatic-guard-cell-homeostasis</footnote>
where $\lambda$ is a Lagrange multiplier of size $N - \dim(\mathcal{Y}_x)$. In (2.26), the first line describes the balance of forces that characterizes the equilibrium point. The multiplier $\lambda$ is adjusted so that the artificial force $W\lambda$ enforces the constraint $\mathbf{n} \in \mathcal{Y}_x$. The second line corresponds to the constraint itself in matrix form.

To obtain a simpler system, we seek to eliminate $\lambda$. The matrix $W$, whose columns form a basis of $\mathcal{Y}_x^\perp$, does not depend on $\mathbf{x}$, as the $\mathcal{Y}_x$ are parallel to each other. We define the matrix $V$, with $N$ rows and $\dim(\mathcal{Y}_x)$ columns, whose columns form a basis of $\mathcal{Y}_x$. In particular, the columns of the block matrix $(V, W)$ form a basis of $\mathbb{R}^N$, and $V^T W = 0$. By applying $V^T$ to the first line of (2.26), we obtain the equivalent system of size $N$

$$ \left\{ \begin{aligned} \mathbf{V}^T \nabla G(\mathbf{n}) &= 0 \\ \mathbf{W}^T (\mathbf{n} - \mathbf{n}^0) &= \mathbf{W}^T \mathbf{A} \mathbf{x}. \end{aligned} \right. \qquad (2.27) $$

This system illustrates well that the coordinates of $\mathbf{n}^*(\mathbf{x})$ along directions collinear to $\mathcal{Y}_x$ are defined by the balance of forces, while the coordinates along directions orthogonal to $\mathcal{Y}_x$ are explicitly determined from $\mathbf{x}$.

Let us write (2.27) in the case of the toy problem from Figure 3, with $V$ and $W$ defined by

$$ \mathbf{V} = \begin{pmatrix} 1 & 0 \\ 2 & 0 \\ 0 & 1 \end{pmatrix} \quad \text{and} \quad \mathbf{W} = \begin{pmatrix} 2 \\ -1 \\ 0 \end{pmatrix}. \qquad (2.28) $$

As the derivatives of $G_{\text{elec}}$ and $G_{\text{mech}}$ do not raise numerical issues, we only specify the expression for the chemical energy derivatives using (2.4) and (2.5), leading to the system

$$ \left\{ \begin{aligned} & \left[ RT \ln \frac{y_{H^+}(\mathbf{n})}{y_{\text{ext,H}^+}} + \frac{\partial(G_{\text{elec}} + G_{\text{mech}})}{\partial n_{H^+}}(\mathbf{n}) \right] + 2 \left[ RT \ln \frac{y_{\text{Cl}^-}(\mathbf{n})}{y_{\text{ext,Cl}^-}} + \frac{\partial(G_{\text{elec}} + G_{\text{mech}})}{\partial n_{\text{Cl}^-}}(\mathbf{n}) \right] = 0 \\ & RT \ln \left( \frac{y_{\text{H}_2\text{O}}(\mathbf{n})}{y_{\text{ext,H}_2\text{O}}} \right) + \frac{\partial(G_{\text{elec}} + G_{\text{mech}})}{\partial n_{\text{H}_2\text{O}}}(\mathbf{n}) = 0 \\ & 2(n_{H^+} - n_{H^+}^0) - (n_{\text{Cl}^-} - n_{\text{Cl}^-}^0) = -2x. \end{aligned} \right. \qquad (2.29) $$

In practice, we aim to solve (2.29) using a Newton method [38, Algorithm 11.1]. Here the presence of the logarithm function in the expression of chemical potentials represents a numerical difficulty, as the logarithm function is only defined for positive parameters, and takes large values for arguments close to zero. This may cause the Newton method to fail, either because the algorithm tries to evaluate the residual somewhere where it is not defined, or because the Jacobian matrix is poorly conditioned. To increase the numerical method robustness, we use a relaxed version of (2.29), following a so-called Cartesian representation technique, recently proposed in [39]. Without entering too much into details, let us mention that, with this technique, the chemical potentials are represented as auxiliary variables $\mu_{H^+}$, $\mu_{\text{Cl}^-}$, $\mu_{\text{H}_2\text{O}}$ instead of functions of $\mathbf{n}$. Then, $N$ new lines are added to (2.29) to enforce a relation of the form
$f(\mu_{i,A}/RT) - g(y_{i,A}(\mathbf{n})/y_{\text{ext,A}}) = 0$, where $f$ and $g$ are functions defined on $\mathbb{R}$ and chosen to exhibit good numerical properties, namely

$$f(u) = \begin{cases} \exp(u) - 1 & \text{if } u \le 0 \\ u & \text{if } u > 0 \end{cases} \quad \text{and} \quad g(v) = \begin{cases} v - 1 & \text{if } v \le 1 \\ \ln(v) & \text{if } v > 1, \end{cases} \qquad (2.30)$$

so that the enforced relation is either

$$\frac{\mu_{i,A}}{RT} = \ln \left( \frac{y_{i,A}(\mathbf{n})}{y_{\text{ext},A}} \right) \quad \text{or} \quad \frac{y_{i,A}(\mathbf{n})}{y_{\text{ext},A}} = \exp \left( \frac{\mu_{i,A}}{RT} \right), \qquad (2.31)$$

depending on the values of $y_{i,A}(\mathbf{n})$ and $\mu_{i,A}$ (note that $y_{i,A}(\mathbf{n})$ may take values out of $[0, 1]$ if $\mathbf{n}$ has negative components). As a consequence, the relaxed system, which is actually solved in the numerical procedure, reads

$$\left\{
\begin{aligned}
& \left[ \mu_{H^+} + \frac{\partial(G_{\text{elec}} + G_{\text{mech}})}{\partial n_{H^+}}(\mathbf{n}) \right] + 2 \left[ \mu_{\text{Cl}^-} + \frac{\partial(G_{\text{elec}} + G_{\text{mech}})}{\partial n_{\text{Cl}^-}}(\mathbf{n}) \right] = 0 \\
& \mu_{\text{H}_2\text{O}} + \frac{\partial(G_{\text{elec}} + G_{\text{mech}})}{\partial n_{\text{H}_2\text{O}}}(\mathbf{n}) = 0 \\
& 2(n_{H^+} - n_{H^+}^0) - (n_{\text{Cl}^-} - n_{\text{Cl}^-}^0) = -2x \\
& f\left(\frac{\mu_{i,H^+}}{RT}\right) - g\left(\frac{y_{i,H^+}(\mathbf{n})}{y_{\text{ext},H^+}}\right) = 0 \qquad (2.32) \\
& f\left(\frac{\mu_{i,\text{Cl}^-}}{RT}\right) - g\left(\frac{y_{i,\text{Cl}^-}(\mathbf{n})}{y_{\text{ext},\text{Cl}^-}}\right) = 0 \\
& f\left(\frac{\mu_{i,\text{H}_2\text{O}}}{RT}\right) - g\left(\frac{y_{i,\text{H}_2\text{O}}(\mathbf{n})}{y_{\text{ext},\text{H}_2\text{O}}}(\mathbf{n})\right) = 0.
\end{aligned}
\right.$$

By solving (2.32) successively for $x$ describing an interval $[0, x_{\text{max}}]$ we can compute the trajectory of the moving equilibrium with a limited risk of numerical instabilities.

In the next section, we apply the quasi-static formalism to perform a simulation involving the two-membrane complex from Figure 1 which symbolizes a guard cell.

## 3 Numerical simulations involving a simplified guard cell model

We now apply the energy-based formalism to a concrete system inspired from plant biology. For this, we consider the ion homeostasis in a guard cell during the stoma opening process. Located on plant leaves and stems, stomata are small pores that allow the plant to exchange water and
carbon dioxide with the atmosphere. Each stoma is surrounded by two guard cells, that control its aperture by changing their shape. When the stoma opens, the change of the guard cell shape is only caused by its volume increase, due to water entering the guard cell. In this section, we simulate the cascade of transport processes, driven by two hydrogen pumps on the plasma and vacuole membranes, which leads to water flowing into the cell. In the beginning of this section, we first describe the parameters that we use to approximate the configuration of a guard cell, along with the choice of transporters. Then, we run a simulation and propose an interpretation of the results.

## 3.1 Description of a minimal model representing a guard cell

The processes leading to stoma opening involve the accumulation in the guard cell of different chemical species in a large amount, including ions and also sugars, to modify the osmotic pressure of the guard cell. The accumulation of ions and sugars depends on a wide variety of transporters localized in cellular membranes [40]. Here, we adapt the proposed strategy to build a minimal model of ion transport in the guard cell, that features a limited number of transporters. To be more specific, we limit the number of transporters and chemical species to those directly involved in the modification of the osmotic load. Therefore, we do not include calcium ions, since their function in guard cells is related to signaling and they do not contribute significantly to the osmotic effects [6, 40]. During stoma opening, water uptake is mostly driven by the accumulation of potassium and anions like chloride and malate in the guard cells [18]. Chloride and malate play similar roles, electrical charge compensation and osmotic load, and therefore we only use chloride as a representative for the anions involved in this process. Last but not least, we include hydrogen and water in the model, as the former is transported by hydrogen pumps that drive the whole process, and the latter is directly responsible for cell volume changes. In summary, the simplified model involves four chemical species: hydrogen ($H^+$), chloride ($Cl^-$), potassium ($K^+$), water ($H_2O$).

The guard cell is represented by the two-membrane complex summarized in Figure 1. It contains two compartments, the cytoplasm (compartment 1) and the vacuole (compartment 2). In this context, the state vector is composed of the eight variables

$$ \mathbf{n} = (n_{1,H^+} \ n_{1,Cl^-} \ n_{1,K^+} \ n_{1,H_2O} \ n_{2,H^+} \ n_{2,Cl^-} \ n_{2,K^+} \ n_{2,H_2O}) \in \mathbb{R}^8, \quad (3.1) $$

where the first four components are the chemical amounts in the cytoplasm and the other four are the chemical amounts in the vacuole.

The model includes a selection of the main transporters involved in the stoma opening process [see for instance 41, Fig. 6], keeping in mind that chloride stands as a general anion here. The system includes one hydrogen pump on each membrane (i.e. the plasma and vacuolar membrane) with the plasma membrane pump ejecting hydrogen from the cytoplasm to the external environment [42] and the vacuolar membrane pump ejecting hydrogen from the cytoplasm into the vacuole [43]. Further, each membrane includes a transporter coupling hydrogen and chloride. In the vacuolar membrane, the hydrogen/chloride transporter is an antiporter, as it transports the two species in opposite ways. The molecular identity of this antiporter is known [44, 45, 46]. In the plasma membrane, we introduce a chloride/hydrogen symporter
<table><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Cell diameter</td><td>10 μm</td></tr><tr><td>Cell wall surface stiffness</td><td>30 Pa·m</td></tr><tr><td>Initial pressure</td><td>0.5 MPa</td></tr><tr><td>Initial plasma membrane potential</td><td>-50 mV</td></tr><tr><td>Initial vacuole membrane potential</td><td>-40 mV</td></tr></tbody></table>

<table><caption>Table 1: Numerical values of chosen parameters and initial conditions for the simulation</caption><thead><tr><td>Parameter</td><td>External</td><td>Cytoplasm</td><td>Vacuole</td></tr></thead><tbody><tr><td>pH</td><td>5.7</td><td>7.1</td><td>5.8</td></tr><tr><td>KCl concentration</td><td>30 mmol/L</td><td>15 mmol/L</td><td>20 mmol/L</td></tr><tr><td>Compartment volume</td><td>—</td><td>785 μm<sup>3</sup></td><td>3142 μm<sup>3</sup></td></tr></tbody></table>

(which transports both species in the same direction), that has not yet been molecularly identified [47, 40]. For this reason, its stoichiometry is not known. To obtain symmetric behaviors between transporters at the plasma and vacuolar membrane, we chose the following coefficients for this symporter: $H^+ + 2Cl^-$. Using different stoichiometric coefficients for the plasma membrane hydrogen/chloride symporter, such as $H^+ + Cl^-$ or $H^+ + 3Cl^-$, changes the amount of osmolytes imported in the cell per pump cycle, but it does not invalidate the analysis below. Finally, both membranes include a potassium channel, which is a passive transporter allowing this ion species to circulate, and water can freely circulate between compartments. For more details on ion transport in guard cells, see the review paper by Roelfsema and Hedrich [18].

An advantage of using a quasi-static model is that it relies on few parameters, as all parameters associated with the transporter dynamics are not required. However, it is important to choose realistic values for the remaining parameters. The parameter values are gathered in Table 1. They mostly revolve around the geometrical properties of the cell model, as well as initial values for the chemical amounts and potentials. First, cell dimensions define the amount of water in each compartment. The guard cell is modeled as a cylinder that elongates in the longitudinal direction, where the diameter and initial length are about 10 µm and 50 µm, respectively [21, Table 1]. Following Mirasole et al. [48, Figure 3], we assume that the vacuole accounts for 30 percent of the whole cell volume as the simulation begins. From these geometric parameters, we can evaluate the initial amount of water in both compartments.

Another parameter attached to the cell geometry is the stiffness of the wall, which is mainly a function of its Young modulus and its thickness. Several studies converge to agree on a value of $E \approx 100$ MPa for the Young modulus in the longitudinal direction and $\delta \approx 0.1 - 1$ µm for the wall thickness [22, 23, 24], which should result in a surface stiffness modulus $E_s = E \times \delta \approx 10 - 100$ Pa·m. In the computational model, we set the surface stiffness modulus to 30 Pa·m, which is consistent with previous studies while giving suitable volume changes in the simulation results.

We choose initial conditions in accordance with the literature. In their review, Jezek and Blatt [40] gather values for ion concentrations and potentials for open or closed stomata. Following their numerical values, we set the initial pressure at 0.5 MPa [see also 49], while the initial
electric potential across the plasma membrane and the tonoplast are set to -50 mV and -40 mV, respectively, using the convention from [50]. Concerning the initial concentrations for hydrogen, potassium and chloride, we set them to realistic values in the context of *in vitro* experiments involving guard cells [41, 48].

## 3.2 Numerical simulations involving the minimal guard cell model

We now show simulations results involving the guard cell model from Figure 1. The system input is composed of the progress of each one of the two hydrogen pumps. These two input variables are stored in the vector $\mathbf{x} = (x_1, x_2) \in \mathbb{R}^2$, with $x_1$ denoting the progress (in mol) of the plasma membrane pump, and $x_2$ denoting the progress of the vacuole membrane pump. As mentioned in Section 2.2, $x_1$ and $x_2$ will both be functions of a scalar variable $s$, in order to obtain an evolution along a one-dimensional trajectory.

The first simulation involves a fixed ratio $\alpha > 0$ between the two pump rates, i.e. $\mathbf{x}(s) = (s, \alpha s)$. For each new value of $s$, the simulation computes a corresponding equilibrium state $\mathbf{n}^*(\mathbf{x}(s))$. Figure 5 shows the evolution of chemical amounts as a function of $s$, along with other physical quantities that can be computed from them. The figure includes pH values and the charge imbalance $Q_i(\mathbf{n})$ across membranes, while the hydrostatic pressure and electric potentials across membranes are defined by

$$P(\mathbf{n}) = \frac{1}{v_{\text{H}_2\text{O}}} \frac{\partial G_{\text{mech}}}{\partial n_{1,\text{H}_2\text{O}}} (\mathbf{n}) \quad \text{and} \quad \phi_i(\mathbf{n}) = \frac{Q_i(\mathbf{n})}{C_i(\mathbf{n})}, \qquad (3.2)$$

respectively (see (2.12) and (2.8) for notations). Electric potentials are plotted following the sign convention from [50], i.e a negative potential value across a membrane means that the cytoplasm is at a lower electric potential than the other compartment. The simulated system generally follows the biologically expected behavior, as the activation of the pumps results in chloride, potassium, and finally water entering the cell. The turgor pressure varies between 0.5 and 4.2 MPa, while the cell volume increases from 3900 to 5100 $\mu\text{m}^3$ (given a molar volume of water of $1.8 \cdot 10^{-5}\text{m}^3/\text{mol}$), in agreement with experimental values reported in [49, Figure 4]. Electric potentials evolve from -50 mV to -150 mV at the plasma membrane and from -40 mV to -47 mV at the vacuole membrane, which is consistent with ranges from [40, Table I]. Finally, concentrations in chloride and potassium in both compartments increase from 15 – 20 mmol/L to ~ 775 mmol/L. These values are in the same order of magnitude as values from [40, Table II], keeping in mind that chloride represents a generic anion species in our model. Given the model simplicity, variables changing with the right amplitude and in the right direction is probably the best we can expect in terms of comparison with experimental data. These ranges of variations will be at the center of our attention in the coming sections. Other aspects of the model behavior do not necessarily reflect the behavior of real guard cells. The behavior of hydrogen concentrations, in contrast, is less expected. While the pH usually remains stable in the cytoplasm and slightly acidic in the vacuole, in our simulations both compartments become strongly basic as the cell gets emptied of its hydrogen. This point is addressed in the next section. However, it should be noted that the hydrogen depletion is very slow. At the end of the simulation, the amount of hydrogen in the cell has decreased by $2 \cdot 10^{-18}\text{mol}$, while the amount of hydrogen actively ejected by the pump is $2 \cdot 10^{-12}\text{mol}$.
<figure><img src="image_6.png" /><figcaption>Figure 5: Simulation results for the guard cell model. The simulation output consists of the evolution of the number of moles of each chemical species as the progress of active transporters increases. Number of moles of water, chloride and potassium are plotted in a cumulated way, with the upper part representing the number of moles in the cytoplasm and the lower part representing the number of moles in the vacuole. While pH reflects the hydrogen concentration, pressure and electric potentials are obtained from derivatives of the energy function.</figcaption></figure>
In this configuration, numerical experiments suggest that the plasma membrane pump drives ion and water transport across the plasma membrane, while the vacuole membrane pump drives exchanges between the cytoplasm and the vacuole. For instance, in the results from Figure 5, the ratio between pump rates is $\alpha = 50\%$, which means that when the plasma membrane pump ejects 100 protons from the cell, the vacuolar pump sends 50 protons into the vacuole. As a result, 50% of the osmolytes entering the cell are sent to the vacuole, which ends up representing 50% of the whole cell volume. To confirm this trend, we run the same simulation with other values for $\alpha$, namely 40%, 80% and 100%. We do not explore values of $\alpha$ greater than 100%, as they raise numerical instabilities due to chemical amounts in the cytoplasm becoming close to zero. Results, gathered in Figure 6, show that, at the end of the simulation, the fraction of solute and water in the vacuole asymptotically matches the ratio between the pump rates. If the stoichiometric coefficients of the hydrogen/chloride symporter in the plasma membrane were changed, then the final amount of osmolytes in the vacuole would not be the same, but it would still be proportional to the vacuolar membrane pump rate. Note that pH, electric potential and pressure values do not change significantly between cases from Figures 5 to 7.

Therefore, the fraction of the cell volume occupied by the vacuole can be controlled through the ratio between the two pump rates. When inducing stoma opening in *in vitro* guard cells, Mirasole et al. [48] observed that the cytoplasm volume remains approximately constant along the process. To reproduce this phenomenon in the simulation, we set as previously the plasma membrane pump rate at 1 (i.e. $x_1(s) = s$), and then we adjust the vacuole pump progress $x_2(s)$ at each simulation step, so that the amount of water in the cytoplasm remains constant. The results of the simulation with constant cytoplasm volume are shown in Figure 7. The vacuole pump rate $x_2'(s)$ can be read as the slope of the vacuole pump progress in Figure 7b. For each value of the progress $s$, it is equal to the current fraction of cell volume occupied by the vacuole (Figure 7a). Initially at 30%, this fraction has increased to 50% by the end of the simulation, while the slope of the vacuole pump progress graph is slightly larger at the end of the simulation than at the beginning.

Because the transport that occurs at one membrane does not influence the transport at the other membrane, it is possible to describe separately the events happening at each membrane, using Figure 1. At the plasma membrane, hydrogen expelled from the cell by the pump is almost completely compensated for by hydrogen entering back through the symporter, accompanied by chloride. This explains why the amount of chloride in the cell is approximately twice as large as the pump progress $s$ (Figure 5c). However, this cycling is not perfect, and the amount of hydrogen entering the cell through the symporter is slightly smaller than that ejected by the pump. Consequently, there is a net loss of hydrogen by the cell. (Figure 5a). The massive intake of chloride is accompanied by an approximately equal intake of potassium (Figure 5d), through the potassium channel. As these two ions represent most of the charge transfer, the difference between chloride and potassium intake can be read in the charge imbalance plot (Figure 5h). Finally, water intake (Figure 5b) is synonymous with volume increase, and the reaction of the elastic cell wall to this deformation can be seen in Figure 5f through the pressure increase. A similar chain of events occurs at the vacuolar membrane, except that hydrogen is pumped into the vacuole, and it goes back into the cytoplasm through the antiporter, in exchange for chloride entering the vacuole. Again, the cycling of protons is not ideal and there is a net loss of hydrogen
<figure><img src="image_7.png" /><figcaption>Figure 6: Simulation results with constant ratios between pump rates. The rates of the vacuole pump are $\alpha = 40\%$ (left), $\alpha = 80\%$ (middle), and $\alpha = 100\%$ (right). The plotted values are the numbers of moles of chloride (top), potassium (middle) and water (bottom), with the amount of water being proportional to the compartment volume. All plots on a same row share the same y-axis. The ratio between osmolyte and water transport at the vacuole membrane and at the plasma membrane is asymptotically similar to the ratio between the hydrogen pump rates, which corresponds to equal concentrations of osmolytes between compartments.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Simulation result involving a constant-volume constraint. The rate of the vacuole membrane pump is controlled to keep the cytoplasm volume constant. The plot (b) shows the progress of both pumps along the simulation. This is an example of tuning the model input parameters to qualitatively reproduce experimentally observed phenomena.</figcaption></figure>
from the vacuole. This, together with the loss of hydrogen at the plasma membrane, results in an increase of the cell pH. An interpretation of the physical processes causing this evolution is proposed in Section 4.3.

### 3.3 Impact of hydrogen buffering mechanisms on simulation results

So far, we have considered a model involving solute and water transport only, without chemical reaction. While simulation results show a plausible evolution for most physical quantities, a discrepancy remains between the simulation and experimentally observed guard cell behavior. In the simulation, hydrogen tends to disappear from the cell, leading to a very basic pH in both compartments. On the other hand, experimental observations suggest that the cytoplasm pH remains close to 7, while the vacuole pH remains acidic or slightly increases, depending on reports [51, 6, 52]. A key parameter in the cell, pH is actually well controlled using a range of mechanisms, called hydrogen buffering. These mechanisms are able to capture, release or transport hydrogen, to maintain the pH in a desired range in each compartment [53]. To remedy this inconsistency in our simulation, we now include a hydrogen buffering mechanism in the model, and we investigate its consequences on the simulation output. Though both the cytoplasm and the vacuole of guard cells feature hydrogen buffering in a guard cell, the vacuole pH is less tightly controlled compared with the cytoplasm pH [51, 6, 52]. For the sake of simplicity, we only consider hydrogen buffering in the cytoplasm.

In reality, hydrogen buffering results from several complex mechanisms, including chemical components (e.g. organic acids, phosphates), metabolic activity, solute transport and buffering capacity of the cell wall [53, 54]. In the model, we summarize these mechanisms in the form of a chemical reaction that compensates the variations of hydrogen concentration in a compartment. We represent buffering as an acid-base reaction of type

$$AH_m \rightleftharpoons A^{m-} + mH^+, \qquad (3.3)$$

where A is a species able to bind with hydrogen. When the hydrogen concentration decreases, the reaction provides some hydrogen by consuming $AH_m$, while the reaction produces $AH_m$ when hydrogen is present in excess. The stoichiometric coefficient m is introduced as a numerical tool to adjust the strength of pH regulation. We choose a large m to create a strong buffering capacity in the cytoplasm, without inducing perturbations of osmotic effects due to the presence of large amounts of buffer molecules in the compartment. Here, we choose m = $10^4$ and we set the initial concentration of the buffering solution to 1 mmol/L in the cytoplasm. In the same fashion as passive transporters, this chemical reaction is represented by one of the vector directions that span $\nabla_x$, as it gives the system the freedom to evolve toward production or consumption of $AH_m$.

Figure 8 shows the simulation results for the buffered model, with the same pump rates as in Figure 5, i.e. 1 for the external pump and 0.5 for the vacuole pump. The hydrogen amounts evolve differently compared with previous simulations. Here, the cytoplasm pH remains stable and close to 7, and the vacuole pH remains in acidic values, which is in agreement with experimental observations. Note that the pH in the vacuole decreases more than what is suggested by experimental measures. The amplitude of pH variations in the vacuole can be corrected
by adding a buffering mechanism in this compartment. Interestingly, the evolution of chloride and potassium amounts is also changed, and the connection between compartment volumes and pump rates observed in Figure 6 is now more loose. While, without buffer, each pump is responsible for making osmolytes cross the membrane it is located on, the introduction of the buffering mechanism makes the effect of both pumps more intricate. Our simulations suggest that the vacuolar membrane pump triggers osmolyte transport all the way from the external environment into the vacuole, and the amount of osmolytes crossing each membrane is not proportional to the pump rates anymore. Though studying the buffered system in details is out of the scope of this article, results from Figure 8 show that our model can be extended to produce more realistic simulations.

# 4 Sensitivity analysis using second derivatives of energy

In this section, we perform a sensitivity analysis of the guard cell model, based on second-order information about the energy function. Our aim here is to explain the range of variation predicted by the model for all variables, and also to highlight the couplings between variables created by the different terms of the energy function. After a short calculation that emphasizes the role of the energy second derivatives in determining the system response to a perturbation, we give some insight about the physical meaning of the energy second derivatives. Finally, we implement this reasoning by comparing second derivative coefficients and simulation results in the case of the guard cell model. As G is a regular function of several variables ($n_1, n_2, \dots$), its second derivatives at $\mathbf{n}$ are stored in a symmetric matrix called the Hessian matrix, denoted by $\nabla^2 G(\mathbf{n})$, and the Hessian coefficient at index $(i, j)$ is denoted by $\frac{\partial^2 G}{\partial n_i \partial n_j}(\mathbf{n})$.

## 4.1 Second derivatives of the energy function: a physical analogy

First, let us clearly identify where the energy second derivatives appear in the calculations. Remember that, for a fixed $\mathbf{x}$, the equilibrium state $\mathbf{n}^*(\mathbf{x})$ is solution to the nonlinear system (2.27). Now, assume that $\mathbf{x}$ is subject to a perturbation $d\mathbf{x}$. The equilibrium perturbation $d\mathbf{n}^*$ is the solution to a tangent linear system involving the Hessian matrix of the energy function. This linear system, obtained by differentiating (2.27), reads

$$ \left\{ \begin{array}{l} \mathbf{V}^T \nabla^2 G(\mathbf{n}) \, d\mathbf{n} = 0 \\ \mathbf{W}^T \, d\mathbf{n} = \mathbf{W}^T \mathbf{A} \, d\mathbf{x}. \end{array} \right. \qquad (4.1) $$

Like previously, the first line of (4.1) concerns the coordinates of $d\mathbf{n}$ in the direction of $\mathcal{V}_x$, which adjust to re-establish a physical equilibrium, while the remaining coordinates are explicitly determined from $d\mathbf{x}$ by the second equation. The Hessian matrix plays a role in the evolution of $\mathbf{n}$ in the directions of passive transporters, as a reaction to an imposed change in the directions orthogonal to $\mathcal{V}_x$. While forces are defined by the first derivatives of the energy, this analysis explains why second derivatives must be considered in a sensitivity analysis as they reflect force responses to perturbations. This idea can be illustrated through a simple physical analogy. Consider a spring of stiffness $k$, clamped at one extremity and subject to a force $F$ at the other
<figure><img src="image_9.png" /><figcaption>Figure 8: Simulation result for the guard cell model with hydrogen buffer. The plotted physical quantities are the same as in Figure 5. The cytoplasm pH values follow a more realistic behavior when hydrogen buffering is taken into account.</figcaption></figure>
extremity, causing it to extend by a length x in a one-dimensional way. The spring elastic energy function reads

$$E(x) = \frac{1}{2}kx^2. \qquad (4.2)$$

Also, $E'$ and $E''$ denote the first and second derivatives of the single-variable function $E$. Here the second derivative of $E$ is the stiffness $k$, which defines the spring resistance to changes of lengths. If the force applied by the external user at the spring extremity varies from $F$ to $F + \Delta F$, the spring length changes by $\Delta x = \Delta F/k$. In other words, the spring stiffness gives an indication about how much the system resists to perturbations. In the same fashion, if we consider a nonlinear spring with an energy function $E$, the force applied at its extremities satisfies $F = E'(x)$, and as a consequence, a force perturbation $\Delta F$ causes the length to change at first order by $\Delta x = \Delta F/E''(x)$, where $E''(x)$ plays the role of an immediate stiffness.

The second derivatives of $G$ are like the stiffness of a spring. However, as opposed to the spring energy function, $G$ is a function of several variables, and its second derivatives are gathered in the Hessian matrix $\nabla^2 G(n)$. In this matrix, nondiagonal terms denote couplings between variables, where moving one variable results in a force causing the variation of another variable. An example of such coupling appears when a transfer of charged particles affects charged particles from another species. In the next section, this difficulty is mitigated by computing the eigenvectors of $\nabla^2 G(n)$, which allows to identify decoupled physical balances in the system.

Pushing the spring analogy a little further, we now consider $n$ springs in series. The spring chain is clamped at both ends, so that its total length is imposed. At equilibrium, the individual deflection of the $i$-th spring is denoted by $x_i$, and due to the imposed total length, the spring deflections all add up to an imposed value $y$. The $x_i$ values are determined by solving the minimization problem

$$\min_{x_1, \dots, x_n} \frac{1}{2} k_1 x_1^2 + \dots + \frac{1}{2} k_n x_n^2 \quad \text{subject to constraint} \quad x_1 + \dots + x_n = y, \qquad (4.3)$$

where $k_i$ is the stiffness of the $i$-th spring. Now, if the imposed cumulated deflection $y$ changes by a perturbation $\Delta y$, the system will find a new equilibrium state, where each deflection $x_i$ has changed by a small displacement

$$\Delta x_i = \frac{1/k_i}{1/k_1 + \dots + 1/k_n} \Delta y. \qquad (4.4)$$

The sensitivity of each spring to the perturbation is inversely proportional to its stiffness. If the first spring has the largest stiffness, then its deviation from its initial state is smaller than the deviation of other springs, but it is also even smaller as there are many springs in the chain to absorb the perturbation. In other words, we could say that the initial state of springs with a larger stiffness is preserved *at the expense* of that for springs with a lower stiffness. Similarly, if each spring is nonlinear and characterized by its energy function $E_i(x_i)$, then, solving a minimization problem similar to (4.3) yields the first-order responses to a perturbation $\Delta y$,

$$\Delta x_i = \frac{1/E_i''(x_i)}{1/E_1''(x_1) + \dots + 1/E_n''(x_n)} \Delta y. \qquad (4.5)$$

The balances that exist in our system, such as the balances of concentrations, of charges, or the mechanical equilibrium of the system, can be compared with the springs in series, with the
constraint that the springs have their extremities attached together corresponding to the transport constraints imposed by the stoichiometry of transporters. Each part of the energy function promotes one balance, taking higher values when the balance is not respected. But, somehow, all these balances are in competition with each other, and, when a perturbation occurs, the term with the largest second-order derivative wins, i.e. the corresponding balance is more preserved, at the expense of other balances. In the remaining of this section, we exploit this analogy to interpret the simulation results from Section 3.

## 4.2 Comparison between Hessian matrices and simulation results

In this section, we compare the numerical properties of the energy Hessian matrix the range of variation of chemical amounts in the simulation results. While this analysis is straightforward for systems involving a few variables, the exercise becomes increasingly difficult when more variables are added, because Hessian matrices become larger. To illustrate the role of Hessian coefficients, we focus on the guard cell model without hydrogen buffering, as it features fewer variables and is easier to interpret than the model with hydrogen buffering.

Figure 9 shows the Hessian matrix for each term of the energy function, along with the Hessian matrix for the total energy function. These matrices are evaluated at the final state $n^*(x_{max})$, but their structures remain similar all along the simulation (see Hessian matrices evaluated at the initial state in S1 Fig.). Concerning the axis labels, each row and column of a Hessian matrix is associated with one coordinate of $n$. Chemical amounts are organized in the same order as in (3.1), i.e. amounts in the cytoplasm first and amounts in the vacuole next. For instance, the coefficient at row 2 and column 7 corresponds to the coupling between $n_{1,Cl}^-$ and $n_{2,K^+}$. The values written in Figure 9 are in logarithmic scale, i.e. 27.2 on the matrix heatmap means that the corresponding coefficient is equal to $10^{27.2}$ J/mol$^2$ in absolute value, while a blank box denotes a null coefficient. Therefore, a number visible in the total Hessian heatmap (Figure 9d) is not the sum of the corresponding numbers in the other heatmaps. Instead, it often represents the largest order of magnitude among the corresponding coefficients from the other matrices.

To investigate the important balances in the system, we plotted in Figure 10a the eigenvectors $u_0(s_{max}), \cdots, u_8(s_{max})$, and eigenvalues $\lambda_1(s_{max}) < \cdots < \lambda_8(s_{max})$ of the total energy Hessian matrix, evaluated at the final state. Figure 10b shows the evolution of eigenvalues during the simulation. Except for $\lambda_3(s) \approx \lambda_4(s)$, the eigenvalues remain well separated and keep the same ordering. In Figure 10c, we plot the angle between each eigenvector $u_i(s)$ and its final direction $u_i(s_{max})$. Since eigenvectors are unit vectors, a small angle between two of them means that they are almost equal to each other. Here, for most of the simulation, eigenvectors stay less than one degree away from their final direction. We thus consider that the Hessian matrix eigenvectors remain approximately constant and equal to the final eigenvectors, with only eigenvalues changing along simulation steps, and we will base our sensitivity analysis on the final eigenvectors.

Working with eigenvectors allows to break the space of states into a set of orthogonal directions $u_i$, each one representing a physical process. The system resists to displacements in each direction $u_i$ independently from the others, to an extent given by the eigenvalue $\lambda_i$. It is striking
<figure><img src="image_10.png" /><figcaption>Figure 9: Hessian matrices for each term of the energy function, evaluated at the final state of the simulation. Hessian coefficients (in J/mol²) are represented in logarithmic scale, i.e. 27.2 means 10²⁷.² and a blank case means that the coefficient is zero. The Hessian matrices define the sensitivity of chemical, electrostatic and mechanical forces to perturbations.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Spectral information about the total energy Hessian matrix. (a) Eigenvectors $\mathbf{u}_1, \dots, \mathbf{u}_8$ of the total energy Hessian at the final state, with corresponding eigenvalues $\lambda_1, \dots, \lambda_8$. (b) Evolution of the Hessian eigenvalues along the simulation steps. Eigenvalues change but keep the same ordering during the process, and the corresponding eigenvectors evolve continuously during the simulation. (c) Angle between each eigenvector $\mathbf{u}_i(s)$ at a given step and its final position. A small angle means that eigenvectors are close to each other. After a short transient regime, eigenvectors remain less than one degree apart from their final position.</figcaption></figure>
that most eigenvectors only contains one or two nonzero coordinates, resulting in independent
processes involving one or two species. By looking at the position of their nonzero coordinates,
eigenvectors can be connected with the physical balances enforced by the energy function. Also,
many eigenvalues appear as coefficients in the chemical, electrostatic or elastic energy Hessian
matrix. Based on this observation, we can connect each eigenvector to the energy term that is
responsible for enforcing the corresponding balance.

First, the two largest eigenvalues $\lambda_7$ and $\lambda_8$ appear in the chemical energy Hessian matrix as diagonal parameters, as well as $\lambda_3$ and $\lambda_4$. The eigenvectors $\mathbf{u}_7$ and $\mathbf{u}_8$, correspond to variations of hydrogen amounts in the vacuole and cytoplasm, respectively, while $\mathbf{u}_3$ and $\mathbf{u}_4$ are connected with equal variations of chloride and potassium amounts in the vacuole and cytoplasm, respectively. In the chemical Hessian matrix, the diagonal coefficients associated with the amounts of ion species are much larger in magnitude than nondiagonal parameters. The coefficient associated with hydrogen is approximately $10^{14}$ times larger than nondiagonal parameters, while the coefficients associated with chloride and potassium are ~ $10^2$ times larger. This suggests that chemical forces mainly act on the individual evolution of each one of these amounts of substance, without creating a significant coupling between ion species. The chemical energy term influences the variations of chloride and potassium independently, but in the total energy function it only has a visible impact on equal variations of these two species. In contrast, opposed variations of chloride and potassium amounts have an impact on the balance of charges, which is mostly influenced by electrostatic forces.

Eigenvectors $\mathbf{u}_5$ and $\mathbf{u}_6$ correspond to physical processes where chloride and potassium amounts evolve in opposite directions, resulting in charge transfers across membranes. Moving in the direction $\mathbf{u}_5$ results in transferring charges from the cytoplasm to the vacuole and external environment, while in $\mathbf{u}_6$, charges are moved from both compartments to the external environment. Also, $\lambda_5$ appears in the electrostatic Hessian matrix, and $\lambda_6$ is close to certain coefficients in the same matrix (see framed zones in Figure 9b). The nondiagonal structure of the electrostatic Hessian matrix comes from the two terms of the electrostatic energy expression (2.8). One term, associated with the plasma membrane, depends on amounts in both compartments, while the other term, associated with the vacuole membrane, only depends on amounts in the vacuole. In the matrix, the three identical 4 × 4 blocks are due to the first term, while the bottom right block, different from the others, includes in addition second derivative information from the second term. Though the framed zones in the matrix indicate that electrostatic forces affect all ion species, forces influencing hydrogen amounts are dominated by chemical effects (see $\mathbf{u}_7$, $\mathbf{u}_8$ in Figure 10a), which explains why the two eigenvectors related to electrostatic effects only involve chloride and potassium.

The two remaining eigenvectors are connected with movements of water between compartments, with different eigenvalues depending on whether volume change is involved or not. The direction $\mathbf{u}_2$ represents same-sign variations of water amounts in both compartments, resulting in variations of the cell volume. The corresponding eigenvalue $\lambda_2$ is close to coefficients appearing in the elastic Hessian matrix. Even though the chemical Hessian matrix exhibits diagonal coefficients for water that are in the same order of magnitude, these coefficients are smaller in the initial Hessian matrix (S1 Fig. a). The fact that $\lambda_2(s)$ remains constant like the elastic Hessian matrix suggests that the term offering the most resistance to variation of cell volume is the elastic
energy term. In this case, the volume of water entering the cell as a result of osmotic pressure can be entirely adjusted in our simulation, just by tuning the cell wall stiffness parameter, with a negligible impact on other variables. The last eigenvector **u**₁ represents the exchange of water between the cytoplasm and the vacuole. This exchange of water is not influenced by mechanical effects since it does not involve cell volume changes. Involved physical effects may include chemical forces, but also electrostatic forces, as variations of the vacuole volume have an impact on its surface area, and thus on its electrical capacitance. Actually, no energy term significantly restrains exchanges of water between compartments, since the corresponding eigenvalue λ₁ is twelve orders of magnitude smaller than λ₂, explaining why the ratio between compartment volumes strictly follows that of osmolyte concentrations in our simulations. Though, a positive λ₁ has the merit to ensure the energy function convexity, and, therefore, the stability of the computed equilibria.

Just as the total energy Hessian in Figure 9d gathers the salient visual features from the three matrices that compose it, its eigenvectors and eigenvalues reflect the main physical balances influenced by each term of the energy function. It is noteworthy that the range of eigenvalues spreads over 27 orders of magnitude, resulting in a clear hierarchy between the physical balances, from the tight balance of hydrogen concentrations to the water balances, more loosely enforced. These differences between eigenvalues have their counterpart in the amplitude of variations of certain physical quantities in Figure 5. In Figure 5 (a, c, d), the annotations indicate the variation range of each chemical amount (in mol) in the cytoplasm (A) and in the vacuole (B). In particular, the variations of hydrogen amounts (A ≈ 10⁻¹⁹ mol, B ≈ 10⁻¹⁸ mol) are several orders of magnitude smaller than variations of chloride and potassium (A ≈ B ≈ 10⁻¹² mol), which can be explained by the hydrogen variations being highly restricted by a stronger feedback force (λ₈, λ₇ ≫ λ₄, λ₃). Comparing the variations of the same physical quantity between the two compartments is less straightforward. The variation range of hydrogen amounts in both compartments is mostly a consequence of initial conditions, without significant impact of the pump rates. On the other hand, we saw in Figure 6 that the ratio between chloride and potassium quantities in both compartments can be fully explained by the ratio between the pump rates.

To study the electrostatic balances ($\mathbf{u}_5, \mathbf{u}_6$) in the system, we plotted in Figure 5h the charge imbalance across each membrane, defined by $Q_i(\mathbf{n})/F$ (see 2.8). It denotes the amount of charges (in mol) stored in each membrane/capacitor (see the circuit analogy in Figure 2). Like previously, the annotations show the charge imbalance variations at the plasma membrane (A ≈ 10⁻¹⁷ mol) and at the vacuole membrane (B ≈ 10⁻¹⁸ mol). Because charge imbalance shares units with chemical amounts, its variations (Figure 5h) can be compared with the amount variations of all species (Figure 5a, c, d). Variations of the charge imbalance at the plasma membrane are slightly larger in amplitude than variations of the amount of hydrogen in the cell, and much smaller than variations of the amount of chloride and potassium in the cell (10⁻¹⁸ mol < 10⁻¹⁷ mol ≪ 10⁻¹² mol). This is consistent with eigenvalues associated with electrostatic effects being smaller than eigenvalues associated with hydrogen, and larger than eigenvalues associated with chloride and potassium (10²⁶.⁵ ≫ 10²⁰.⁷ ≫ 10¹⁵.¹). Concerning charge imbalance at the vacuolar membrane (B), simulations show that it is strongly impacted by a combination of initial conditions and pump rates, and it cannot be explained using Hessian eigenvalues only.

The connection between the Hessian matrix eigenvalues and the amplitude of variation of
<figure><img src="image_12.png" /><figcaption>Figure 11: Summary of the transport cascade at the plasma membrane in the guard cell model without buffering, showing the hierarchy of forces at work upon a charge imbalance. Priorities in equilibria is decreasing from 1 to 4. The vacuole is not represented in the figure. In steps 2 and 3, one balance is restored at the expense of another less strictly enforced balance.</figcaption></figure>

the system variables illustrates well the role played by the second derivatives of the energy function. Like springs with a large stiffness, balances associated with large eigenvalues (hydrogen, charges) tend to resist to a perturbation, while balances associated with small eigenvalues (chloride, potassium, water) absorb most of it. However, even springs with a very large stiffness deflect a little when subject to a perturbation, and so does the amount of hydrogen, resulting in the hydrogen loss observed in Figure 5a.

Though the system evolution is also influenced by other parameters such as pump rates and transport stoichiometry, studying the Hessian matrix of the energy function allows us to explain trends in the range of variation of physical quantities. These trends only depend on the energy function and the current system state, and they do not depend on the nature of transporters. As a consequence, this analysis remains valid for similar systems involving embedded compartments and osmolyte transport through membranes, including systems where the disposition of transporters is not known with certainty.

## 4.3 An interpretation for turgor pressure building in plant cells

Comparing the sensitivities of all the balances that govern the system allows to qualitatively understand the hierarchy between the physical phenomena at play, and explain why a property is enforced at the expense of another one. In particular, we can interpret the guard cell scenario from Section 3 through the prism of this hierarchy between the terms of the energy function, go-
ing from the activation of hydrogen pumps all the way to water entering the cell. Our simulation results in Section 3.2 suggest that transport processes on both membranes function independently of each other, each pump driving all transport across the membrane it belongs to. For this reason, we show in Figure 11 an interpretation for solute transport at the plasma membrane, while a similar summary can be made for the vacuole membrane transport. In this cascade of transport reactions, each step restores one balance by perturbing another one that is lower in the hierarchy. For instance, in Figure 11, expelling hydrogen from the cell results in hydrogen entering back through the symporter, accompanied by chloride, as deflecting from the chloride chemical balance is much cheaper than deflecting from the hydrogen chemical balance. In the same fashion, the amount of potassium is forced to stick to that of chloride to restore the balance of charges, as this balance is much more sensitive than that of potassium. Membrane polarization, characterized by an electric potential building across a membrane, reflects the deflection of the system from electric neutrality. Finally, the resulting imbalance in chloride and potassium is compensated for by water entering the cell. The final volume change is a compromise between chemical forces and mechanical forces, and it mostly depends on the mechanical parameters of the cell wall. Note that, if potassium could not enter the cell to ensure an electrically neutral flow of osmolytes, a small amount of chloride would still enter the cell at the expense of the charge balance, as the balance of hydrogen is more sensitive than the balance of charges.

The difference of sensitivities illustrates well how, by acting on hydrogen, which is present in very small amount, the cell manages to transfer into its vacuole an amount of water that is several orders of magnitude larger than that of hydrogen. Generating a perturbation on a strictly enforced balance creates a leverage effect. As it is diverted toward other balances, the perturbation gets amplified and involves larger-scale transport, which remains cheaper than maintaining a hydrogen imbalance. Even though transporters play a key role in determining the possible directions of evolution for the system state, the hierarchy between all forces, as well as involved orders of magnitude, are mostly a matter of energy function.

# 5 Discussion and Conclusion

**Summary** In an attempt to simplify the modeling of multi-physics couplings in multi-compartment plant cell physiology, we proposed in this paper a model based on an energy function integrating electrochemical and hydro-mechanical effects. We applied this energetic approach to model solute and water exchanges in a guard cell during stoma opening. By using a quasi-static scenario, we investigated the system properties without relying on transporter dynamics, thus reducing the number of parameters and the model complexity. We performed numerical simulations on the guard cell model with and without hydrogen buffering mechanisms, and we showed that the numerical values returned by the model are in agreement with the literature. Thanks to an analysis of the system response to perturbations, we were able to reveal the role and influence of each physical force in the system behavior, which is the central contribution of our work.

**Energy-based formulation** In plant cell homeostasis as in many multi-physics systems, using an energy function is a good way to combine the effects of various physical phenomena. Once
the complexity of the system physics is embedded in the energy function expression, the result-
ing system has a clear structure, and the role of each physical effect or transporter is identified.
This clear structure is also reflected when it comes to implementing this variational formulation
on a computer. In particular, using modern automatic differentiation software, all derivatives of
the energy function are automatically evaluated under the hood.

The energy approach makes a distinction between model-specific and generic features, which both contribute to the system evolution. Transport rules, defined by the configuration of transporters, are specific to the proposed plant cell model, but the energy function is more generic. As the sensitivity analysis is based on the energy function only, a similar hierarchy between forces can be expected in other plant cell models. In contrast, formulations using directly ordinary differential equations [28, 9, 6] do not explicitly separate these two contributions.

**Interpretability of the model** One of the main aims of our modeling strategy was to promote the mechanistic interpretability of the model, in order to obtain insight into the functioning of living systems. We chose to build a system that acts as a mathematical function, which takes the progress of active transporters as input variables, and returns an output equilibrium state. Thus, we were able to control the rate of the hydrogen pumps present in each membrane and investigate the effect of each active transporter on the model evolution. Interestingly, in the model the relative activities of the plasma membrane and vacuolar membrane hydrogen pumps were the main parameters defining the relative volumes of the vacuolar and cytosolic compartments in the cell. Further, it was possible to define a coordination strategy between the pumps in both membranes to keep the cytoplasmic volume constant while the turgor pressure increases in the cell, as experimentally observed in guard cells during stoma opening [48]. In summary, the model indicates that the activities of the plasma membrane and vacuolar membrane pumps need coordination to achieve biological processes involving changes of volume of the subcellular compartments.

Following our strategy, we performed a sensitivity analysis of the various forces (electrical, chemical and mechanical) acting in our model, and we ranked these forces based on their influence on the system behavior. This gave relevant insight on the functioning of cells with nested membranes and strong turgor pressure, like plant cells. The ranking shows that hydrogen-related forces overcome, by orders of magnitude, the other chemical, electrical and hydraulic forces. Notably, in the model hydraulic forces are subordinate to all the other forces. The importance of hydrogen-related forces depends on the orders-of-magnitude lower hydrogen amounts in both compartments compared to the other chemical species. Under these conditions, small changes in the hydrogen amounts (as those generated by the hydrogen pumps) have amplified consequences on the amounts of the other species transported in the systems. This analysis explains why a transport network based on hydrogen pumps is an efficient strategy to create a large osmotic and turgor pressure in plant cells. The active transport of hydrogens by the pumps creates a leverage effect on the transport of other ions, enabling the building of a strong osmotic pressure. Interestingly, land plants, which are among the organisms with the largest osmotic pressure difference with the extracellular media, evolved a transport network based on hydrogen pumps and hydrogen-coupled exchangers and symporters. In contrast, animal cells maintain a low osmotic pressure difference with the extracellular media and do not exploit hydrogen-coupled transport
networks in the plasma membrane, but potassium- and sodium-coupled networks. It is therefore
tempting to speculate that this leverage effect is a reason why land plants evolved a transport
network based on hydrogen, which is present in low amounts in the cellular and extracellular
media.

**Possible extensions of the model** Several extensions can be brought to the proposed model without modifying its principles, as it was the case when adding the buffering mechanism. In this paper, we explored the connections between hydrogen pump activity and the building of a turgor pressure. We considered a cascade of events starting with hydrogen active transport and ending with cell deformation in one dimension. However, this cascade can be extended at both ends. On the upstream side, mechanisms causing the transport of hydrogen may be included in the model. In that case, hydrogen pumps are not controlled by the user, and adenosine triphosphate (ATP), the fuel of active hydrogen transport, is introduced as a new species. On the downstream side, more realistic geometric deformation models could be used [see for instance 22], leading to more complex terms for the mechanical energy.

**Quasi-static limit and possible generalization to dynamical model** We chose a quasi-static approach to model the system evolution in a simple way, at the expense of a fine capture of the transport dynamics. This approach is convenient for models with few transporters, but other choices should be considered to describe more complex transporter configurations. Strong mathematical hypotheses should indeed be satisfied so that the output equilibrium state does not remain constant whatever the value of input variables. In the particular case of Figure 4, a necessary condition for $n^*(x)$ to change as a function of $x$ is that the subspaces $\mathcal{V}_x$ and $\mathcal{V}'_x$ be distinct for two values $x \neq x'$, i.e. the passive transporters do not allow the system to come back to its initial state $n^0$. In contrast, if the direction $v_1$ is collinear to $a$ in the figure, then $n^0$ is an element of $\mathcal{V}_x$ for every $x$ in $[0, x_{max}]$. Therefore, the energy minimizer in $\mathcal{V}_x$ always remains $n^*(x) = n^0$, and the active transporters have no effect on the system. In general, we noticed that a sufficient condition for active transporters to induce an evolution of the equilibrium state is the following: at least one active direction $a$ is not a linear combination of the passive directions ($v_1, v_2, \dots$).

Issues may typically appear when the system involves more transporter than species, with
passive transporters allowing to reverse the effect of active transporters. If active transporters
fail to change the system equilibrium, i.e. the system behavior is entirely due to transporter
dynamics, then a time-dependent model should be considered. A simple way to achieve a dy-
namic model using the energetic approach is given in [9]. With our notation, the flux $J$ across a
passive transporter associated with the direction $v$ is proportional to the energy derivative along
$v$, i.e. $J = -\alpha \nabla G(\mathbf{n})^T v$, where $\alpha$ is a conductance parameter specific to the considered trans-
porter. Compared with a quasi-static model, a dynamic model may be slightly more difficult
to interpret, as its output combines influences from the energy function and transporter conduc-
tance. However, it can be used to create more complex behaviors, for instance by considering
transporters whose dynamics is affected by voltage or pH.
# 6 Supporting material

S1 Fig. Hessian matrices at the initial state of the simulation. This figure is similar to Figure 9, but the Hessian matrices are evaluated at the initial state of the simulation $n^0$. In particular, the analysis of Hessian matrices from Figure 9 still apply to these matrices, as similar differences between coefficient orders of magnitude are present here.

## References

1. Enrico Martinoia, Stefan Meyer, Alexis De Angeli, and Réka Nagy. Vacuolar transporters in their physiological context. *Annual Review of Plant Biology*, 63(Volume 63, 2012):183–213, 2012. ISSN 1545-2123. doi: 10.1146/annurev-arplant-042811-105608.

2. Falco Krüger and Karin Schumacher. Pumping up the volume – vacuole biogenesis in Arabidopsis thaliana. *Seminars in Cell & Developmental Biology*, 80:106–112, 2018. ISSN 1084-9521. doi: 10.1016/j.semcdb.2017.07.008.

3. Sabrina Kaiser and David Scheuring. To lead or to follow: contribution of the plant vacuole to cell growth. *Frontiers in Plant Science*, 11:553, 2020. doi: 10.3389/fpls.2020.00553.

4. Kai Dünser, Maria Schöller, Ann-Kathrin Rößling, Christian Löfke, Nannan Xiao, Barbora Pařázková, Stanislav Melnik, Marta Rodriguez-Franco, Eva Stöger, Ondřej Novák, and Jürgen Kleine-Vehn. Endocytic trafficking promotes vacuolar enlargements for fast cell expansion rates in plants. *eLife*, 11:e75945, Jun 2022. ISSN 2050-084X. doi: 10.7554/eLife.75945.

5. Victoria Vitali, Moira Sutka, Gabriela Amodeo, Osvaldo Chara, and Marcelo Ozu. The water to solute permeability ratio governs the osmotic volume dynamics in beetroot vacuoles. *Frontiers in Plant Science*, 7, 2016. ISSN 1664-462X. doi: 10.3389/fpls.2016.01388.

6. Kunkun Li, Christina Grauschopf, Rainer Hedrich, Ingo Dreyer, and Kai R. Konrad. $K^+$ and pH homeostasis in plant cells is controlled by a synchronized $K^+/H^+$ antiport at the plasma and vacuolar membrane. *New Phytologist*, 241(4):1525–1542, 2024. doi: 10.1111/nph.19436.

7. Matthias Kahm, Clara Navarrete, Vicent Llopis-Torregrosa, Rito Herrera, Lina Barreto, Lynne Yenush, Joaquin Ariño, Jose Ramos, and Maik Kschischo. Potassium starvation in yeast: Mechanisms of homeostasis revealed by mathematical modeling. *PLOS Computational Biology*, 8(6): 1–11, 06 2012. doi: 10.1371/journal.pcbi.1002548.

8. Rosario Astaburuaga, Orlando Daniel Quintanar Haro, Tobias Stauber, and Angela Relógio. A mathematical model of lysosomal ion homeostasis points to differential effects of $Cl^-$ transport in $Ca^{2+}$ dynamics. *Cells*, 8(10), 2019. ISSN 2073-4409. doi: 10.3390/cells8101263.

9. Susanne Gerber, Martina Fröhlich, Hella Lichtenberg-Fraté, Sergey Shabala, Lana Shabala, and Edda Klipp. A thermodynamic model of monovalent cation homeostasis in the yeast saccharomyces cerevisiae. *PLOS Computational Biology*, 12(1):1–23, 01 2016. doi: 10.1371/journal.pcbi.1004703.

10. A. J. Staverman. The theory of measurement of osmotic pressure. *Recueil des Travaux Chimiques des Pays-Bas*, 70(4):344–352, 1951. doi: 10.1002/recl.19510700409.

11. A. J. Staverman. Non-equilibrium thermodynamics of membrane processes. *Trans. Faraday Soc.*, 48:176–185, 1952. doi: 10.1039/TF9524800176.
12. O. Kedem and A. Katchalsky. Thermodynamic analysis of the permeability of biological membranes to non-electrolytes. *Biochimica et Biophysica Acta*, 27:229-246, 1958. ISSN 0006-3002. doi: 10.1016/0006-3002(58)90330-5.

13. J.K.E. Ortega and S.W.J. Welch. Mathematical models for expansive growth of cells with walls. *Math. Model. Nat. Phenom.*, 8(4):35-61, 2013. doi: 10.1051/mmnp/20138404.

14. Euan T Smithers, Jingxi Luo, and Rosemary J Dyson. Mathematical principles and models of plant growth mechanics: from cell wall dynamics to tissue morphogenesis. *Journal of Experimental Botany*, 70(14):3587-3600, 05 2019. ISSN 0022-0957. doi: 10.1093/jxb/erz253.

15. R. Huang, A.A. Becker, and I.A. Jones. Modelling cell wall growth using a fibre-reinforced hypere-lastic-viscoplastic constitutive law. *Journal of the Mechanics and Physics of Solids*, 60(4):750-783, 2012. ISSN 0022-5096. doi: 10.1016/j.jmps.2011.12.003.

16. R. G. M. van der Sman. Hyperelastic models for hydration of cellular tissue. *Soft Matter*, 11: 7579-7591, 2015. doi: 10.1039/C5SM01032B.

17. Euan T. Smithers, Jingxi Luo, and Rosemary J. Dyson. A continuum mechanics model of the plant cell wall reveals interplay between enzyme action and cell wall structure. *The European Physical Journal E*, 47(1):1, Jan 2024. ISSN 1292-895X. doi: 10.1140/epje/s10189-023-00396-2.

18. M. Rob G. Roelfsema and Rainer Hedrich. In the light of stomatal opening: new insights into 'the Watergate'. *New Phytologist*, 167(3):665-691, 2005. doi: 10.1111/j.1469-8137.2005.01460.x.

19. Cornelia Eisenach and Alexis De Angeli. Ion transport at the vacuole during stomatal movements. *Plant Physiology*, 174(2):520-530, 04 2017. ISSN 0032-0889. doi: 10.1104/pp.17.00130.

20. Hugh C. Woolfenden, Alice L. Baillie, Julie E. Gray, Jamie K. Hobbs, Richard J. Morris, and Andrew J. Fleming. Models and mechanisms of stomatal mechanics. *Trends in Plant Science*, 23(9): 822-832, 2018. ISSN 1360-1385. doi: 10.1016/j.tplants.2018.06.003.

21. Tobias Meckel, Lars Gall, Stefan Semrau, Ulrike Homann, and Gerhard Thiel. Guard cells elongate: Relationship of volume and surface area during stomatal movement. *Biophysical Journal*, 92(3): 1072-1080, 2007. ISSN 0006-3495. doi: 10.1529/biophysj.106.092734.

22. Hugh C. Woolfenden, Gildas Bourdais, Michaela Kopischke, Eva Miedes, Antonio Molina, Silke Robatzek, and Richard J. Morris. A computational approach for inferring the cell wall properties that govern guard cell dynamics. *The Plant Journal*, 92(1):5-18, 2017. doi: 10.1111/tpj.13640.

23. Ziv Marom, Ilana Shtein, and Benny Bar-On. Stomatal opening: The role of cell-wall mechanical anisotropy and its analytical relations to the bio-composite characteristics. *Frontiers in Plant Science*, 8, 2017. ISSN 1664-462X. doi: 10.3389/fpls.2017.02061.

24. Hojae Yi, Yue Rui, Baris Kandemir, James Z. Wang, Charles T. Anderson, and Virendra M. Puri. Mechanical effects of cellulose, xyloglucan, and pectins on stomatal guard cells of Arabidopsis thaliana. *Frontiers in Plant Science*, 9, 2018. ISSN 1664-462X. doi: 10.3389/fpls.2018.01566.

25. Ross Carter, Hugh Woolfenden, Alice Baillie, Sam Amsbury, Sarah Carroll, Eleanor Healicon, Spyros Sovatzoglou, Sioban Braybrook, Julie E. Gray, Jamie Hobbs, Richard J. Morris, and Andrew J. Fleming. Stomatal opening involves polar, not radial, stiffening of guard cells. *Current Biology*, 27 (19):2974-2983.e2, 2017. doi: 10.1016/j.cub.2017.08.006.
26. Leila Jaafar, Yintong Chen, Sedighe Keynia, Joseph A. Turner, and Charles T. Anderson. Young guard cells function dynamically despite low mechanical anisotropy but gain efficiency during stomatal maturation in Arabidopsis thaliana. *The Plant Journal*, 118(6):1719–1731, 2024. doi: 10.1111/tpj.16756.

27. Paloma Cubero-Font and Alexis De Angeli. Connecting vacuolar and plasma membrane transport networks. *New Phytologist*, 229(2):755–762, 2021. doi: 10.1111/nph.16983.

28. Adrian Hills, Zhong-Hua Chen, Anna Amtmann, Michael R. Blatt, and Virgilio L. Lew. OnGuard, a computational platform for quantitative kinetic modeling of guard cell physiology. *Plant Physiology*, 159(3):1026–1042, 05 2012. ISSN 0032-0889. doi: 10.1104/pp.112.197244.

29. Zhong-Hua Chen, Adrian Hills, Ulrike Bätz, Anna Amtmann, Virgilio L. Lew, and Michael R. Blatt. Systems dynamic modeling of the stomatal guard cell predicts emergent behaviors in transport, signaling, and volume control. *Plant Physiology*, 159(3):1235–1251, 05 2012. ISSN 0032-0889. doi: 10.1104/pp.112.197350.

30. Yizhou Wang, Maria Papanatsiou, Cornelia Eisenach, Rucha Karnik, Mary Williams, Adrian Hills, Virgilio L. Lew, and Michael R. Blatt. Systems dynamic modeling of a guard cell Cl⁻ channel mutant uncovers an emergent homeostatic network regulating stomatal transpiration. *Plant Physiology*, 160 (4):1956–1967, 10 2012. ISSN 0032-0889. doi: 10.1104/pp.112.207704.

31. Carla Minguet-Parramona, Yizhou Wang, Adrian Hills, Silvere Vialet-Chabrand, Howard Griffiths, Simon Rogers, Tracy Lawson, Virgilio L. Lew, and Michael R. Blatt. An optimal frequency in Ca²⁺ oscillations for stomatal closure is an emergent property of ion transport in guard cells. *Plant Physiology*, 170(1):33–42, 12 2015. ISSN 0032-0889. doi: 10.1104/pp.15.01607.

32. Silvere Vialet-Chabrand, Adrian Hills, Yizhou Wang, Howard Griffiths, Virgilio L. Lew, Tracy Lawson, Michael R. Blatt, and Simon Rogers. Global sensitivity analysis of OnGuard models identifies key hubs for transport interaction in stomatal dynamics. *Plant Physiology*, 174(2):680–688, 04 2017. ISSN 0032-0889. doi: 10.1104/pp.17.00170.

33. P Nelson. *Biological physics*. Freeman, New York, 2004. ISBN 978-0-7167-9897-2.

34. Robert G Mortimer. *Physical chemistry*. Academic Press, third edition, May 2008. ISBN 978-0-080-87859-1.

35. James Keener and James Sneyd. *Mathematical Physiology: I: Cellular Physiology*. Springer New York, New York, NY, 2009. ISBN 978-0-387-75847-3. doi: 10.1007/978-0-387-75847-3.

36. Joseph K. E. Ortega. Augmented growth equation for cell wall expansion 1. *Plant Physiology*, 79 (1):318–320, 09 1985. ISSN 0032-0889. doi: 10.1104/pp.79.1.318.

37. James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL http://github.com/jax-ml/jax.

38. Jorge Nocedal and Stephen J. Wright. *Numerical optimization*. Springer Series in Operations Research and Financial Engineering. Springer, New York, second edition, 2006. ISBN 978-0387-30303-1; 0-387-30303-0.
39. Maxime Jonval, Ibtihel Ben Gharbia, Clément Cancès, Thibault Faney, and Quang-Huy Tran. Parametrization and cartesian representation techniques for robust resolution of chemical equilibria. *Journal of Computational Physics*, 522:113596, 2025. ISSN 0021-9991. doi: 10.1016/j.jcp.2024.113596.

40. Mareike Jezek and Michael R. Blatt. The membrane transport system of the guard cell and its integration for stomatal dynamics. *Plant Physiology*, 174(2):487-519, 04 2017. ISSN 0032-0889. doi: 10.1104/pp.16.01949.

41. Elsa Demes, Laetitia Besse, Paloma Cubero-Font, Béatrice Satiat-Jeunemaitre, Sébastien Thomine, and Alexis De Angeli. Dynamic measurement of cytosolic ph and NO₃ uncovers the role of the vacuolar transporter AtCLCa in cytosolic pH homeostasis. *Proceedings of the National Academy of Sciences*, 117(26):15343-15353, 2020. doi: 10.1073/pnas.2007580117.

42. Sylvain Merlot, Nathalie Leonhardt, Francesca Fenzi, Christiane Valon, Miguel Costa, Laurie Piette, Alain Vavasseur, Bernard Genty, Karine Boivin, Axel Müller, Jérôme Giraudat, and Jeffrey Leung. Constitutive activation of a plasma membrane H⁺-ATPase prevents abscisic acid-mediated stomatal closure. *The EMBO Journal*, 26(13):3216-3226, 2007. doi: 10.1038/sj.emboj.7601750.

43. Melanie Krebs, Diana Beyhl, Esther Görlich, Khaled A. S. Al-Rasheid, Irene Marten, York-Dieter Stierhof, Rainer Hedrich, and Karin Schumacher. Arabidopsis V-ATPase activity at the tonoplast is required for efficient nutrient storage but not for sodium accumulation. *Proceedings of the National Academy of Sciences*, 107(7):3251-3256, 2010. doi: 10.1073/pnas.0913035107.

44. A. De Angeli, D. Monachello, G. Ephritikhine, J. M. Frachisse, S. Thomine, F. Gambale, and H. Barbier-Brygoo. The nitrate/proton antiporter AtCLCa mediates nitrate accumulation in plant vacuoles. *Nature*, 442(7105):939-942, Aug 2006. ISSN 1476-4687. doi: 10.1038/nature05013.

45. Mathieu Jossier, Laëtitia Kroniewicz, Fabien Dalmas, D. Le Thiec, Geneviève Ephritikhine, Sébastien Thomine, Hélène Barbier-Brygoo, Alain Vavasseur, Sophie Filleur, and Nahalie Leonhardt. The arabidopsis vacuolar anion transporter, AtCLCc, is involved in the regulation of stomatal movements and contributes to salt tolerance. *The Plant Journal*, 64(4):563-576, 2010. doi: 10.1111/j.1365-313X.2010.04352.x.

46. Stefanie Wege, Alexis De Angeli, Marie-Jo Droillard, Laetitia Kroniewicz, Sylvain Merlot, David Cornu, Franco Gambale, Enrico Martinoia, Hélène Barbier-Brygoo, Sébastien Thomine, Nathalie Leonhardt, and Sophie Filleur. Phosphorylation of the vacuolar anion exchanger atclca is required for the stomatal response to abscisic acid. *Science Signaling*, 7(333):ra65-ra65, 2014. doi: 10.1126/scisignal.2005140.

47. Fang-Qing Guo, Jared Young, and Nigel M. Crawford. The nitrate transporter AtNRT1.1 (CHL1) functions in stomatal opening and contributes to drought susceptibility in arabidopsis. *The Plant Cell*, 15(1):107-117, 01 2003. ISSN 1040-4651. doi: 10.1105/tpc.006312.

48. Filippo Maria Mirasole, Sara Paola Nastasi, Paloma Cubero-Font, and Alexis De Angeli. Vacuolar control of stomatal opening revealed by 3D imaging of the guard cells. *Scientific Reports*, 13(1): 7647, May 2023. ISSN 2045-2322. doi: 10.1038/s41598-023-34273-x.

49. Peter J. Franks, Thomas N. Buckley, Joseph C. Shope, and Keith A. Mott. Guard cell volume and pressure measured concurrently by confocal microscopy and the cell pressure probe1. *Plant Physiology*, 125(4):1577-1584, 04 2001. ISSN 0032-0889. doi: 10.1104/pp.125.4.1577.
50. Adam Bertl, Eduardo Blumwald, Roberto Coronado, Robert Eisenberg, Goeff Findlay, Dietrich Gradmann, Bertil Hille, Kurt Kohler, Hans-Albert Kolb, Enid MacRobbie, Gerhard Meissner, Christopher Miller, Erwin Neher, Philip Palade, Omar Pantoja, Dale Sanders, Julian Schroeder, Clifford Slayman, Roger Spanswick, Alan Walker, and Alan Williams. Electrical measurements on endomembranes. *Science*, 258(5084):873–874, 1992. doi: 10.1126/science.1439795.

51. Zaida Andrés, Javier Pérez-Hormaeche, Eduardo O. Leidi, Kathrin Schlücking, Leonie Steinhorst, Deirdre H. McLachlan, Karin Schumacher, Alistair M. Hetherington, Jörg Kudla, Beatriz Cubero, and José M. Pardo. Control of vacuolar dynamics and regulation of stomatal aperture by tonoplast potassium uptake. *Proceedings of the National Academy of Sciences*, 111(17):E1806–E1814, 2014. doi: 10.1073/pnas.1320421111.

52. Seoyeon Cha, Wang Ki Min, and Hak Soo Seo. Arabidopsis COP1 guides stomatal response in guard cells through pH regulation. *Communications Biology*, 7(1):150, Feb 2024. ISSN 2399-3642. doi: 10.1038/s42003-024-05847-w.

53. A Kurkdjian and J Guern. Intracellular pH: Measurement and importance in cell activity. *Annual Review of Plant Biology*, 40(Volume 40, 1989):271–303, 1989. ISSN 1545-2123. doi: 10.1146/annurev.pp.40.060189.001415.

54. Jaroslaw Poznanski, Pawel Szczesny, Katarzyna Ruszczyńska, Piotr Zielenkiewicz, and Leszek Paczek. Proteins contribute insignificantly to the intrinsic buffering capacity of yeast cytoplasm. *Biochemical and Biophysical Research Communications*, 430(2):741–744, 2013. ISSN 0006-291X. doi: 10.1016/j.bbrc.2012.11.079.