Multi-physics modeling for ion homeostasis in multi-compartment plant cells using an energy function
Guillaume Mestdagh, Alexis De Angeli, Christophe Godin
TL;DR
This paper introduces a new energy-based model to study how plant cells regulate water and ion transport through multiple compartments.
Contribution
The paper presents a novel energy-based approach to unify chemical, electrical, and mechanical processes in multi-compartment plant cell modeling.
Findings
The energy-based model systematically derives equations for ion and water transport in plant cells.
The model explains the direction of system changes in response to perturbations during stoma opening.
Numerical simulations reveal the role of hydrogen pumps in the transport process.
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…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Fig 1
Fig 2
Fig 3
Fig 4
Fig 5
Fig 6
Fig 7
Fig 8
Fig 9
Fig 10
Fig 11- —Agence nationale pour la recherche
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsPlant nutrient uptake and metabolism · Lipid Membrane Structure and Behavior · Plant and Biological Electrophysiology Studies
1. Introduction
To achieve key living functions, plant cells change their size and shape by regulating the amount of water entering their vacuole [1–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 electrochemical 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 multi-membrane 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–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 [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 [6,27] 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–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 one 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. Methods and models
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 Sect 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 Fig 1, involves two membranes, four chemical species (hydrogen, chloride, potassium, water) and eight transporters.
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).
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 ni,A, and all chemical amounts are gathered in the state vector . 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 Fig 1, of size N = 8, reads
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
where G is an energy function to minimize, is the variable that describes the system state, and is the range of states attainable by the system, which depends on an input parameter . 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 , an electric potential in and a pressure in , but the potential energy is expressed in .
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 the volume enclosed in the i-th membrane, and by the area of its lateral surface. For a cylinder with fixed radius r, the lateral surface area and the volume satisfy . 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 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
where the mole fraction of species A in compartment i, and and 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 ni,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
Here, the external environment is considered as a reservoir, where each species A is present with the constant mole fraction , so that the chemical potential of A in the external environment is constant. The expression for chemical energy is chosen so that , which yields
Note that only depends on chemical amounts in compartment i, so that (5) is actually made of independent terms, with each term depending on mole fractions in one compartment. From expressions (5) and (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 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^2^ [35, Section 2.6].
The electric structure of a multi-membrane complex can be represented by an equivalent circuit, shown in Fig 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 into compartment 1 and into compartment 2, then the charge stored in capacitor C2 changes by , while the charge stored in capacitor C1 changes by + (here, charge displacements act like current intensities). In other words, the charge stored in the capacitor Ci 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 , which does not change the charge stored in capacitor C1.
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 C1 and C2. This circuit analogy is used to evaluate the electrostatic energy term.
Let us now derive an expression for the electrostatic energy stored in capacitor Ci. If we denote by the vector of all chemical amounts in compartment k, the charge contained in compartment k is defined by
where stores the charge of each species, and F is the Faraday constant. The parameters and are fixed as part of the initial conditions. The expression for the charge stored in capacitor Ci reads
where denotes the stored charge in initial conditions. The energy stored in capacitor Ci is defined using the classic formula
Note that the initial membrane charge can be deduced from experimental measures of initial membrane potentials. The expression for the capacitance 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
where is the capacitance per unit surface.
Finally, we evaluate the total electrostatic energy by adding up the energy stored in each capacitor, i.e.
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
where is the molar volume of water. Keep in mind that 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 , 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 , by looking at (11) we obtain
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 . We denote by the initial cylinder lateral surface and by its area. Using a linear constitutive law, the energy integrated over the whole lateral surface reads
where is a surface elastic modulus (in ).
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 [20,22].
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 , which can be a state measured at the beginning of an experiment. However, in a quasi-static model, 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 represents an equilibrium state, we add a linear correction term to the energy function G, where is adjusted so that , which means that is a minimum of G. Namely,
The force 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 . For a given , the total energy function that is actually minimized during the simulation is defined by
Though this energy function depends on , we keep calling it G in the remaining of this section, as is a fixed parameter of the model.
2.2. Membranes, transporters, directions
The energy function described in Sect 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 , 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 in (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 Fig 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.
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.
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 , so that the state vector now reads
The vector 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 Fig 3, the first transporter is a hydrogen pump. It ejects hydrogen from the compartment to the external environment, which results in the hydrogen amount decreasing. As is the first coordinate of the state vector , this is represented in the state space by evolving along the vector
where we have used the same ordering as (16) for chemical species. Once the hydrogen pump has ejected x moles of hydrogen (and no other transport occurred), the new state now reads
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 and , are defined after the stoichiometry of the two transporters (see Fig 3), namely
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 must remain constant. Here, the state can only move along two directions and , 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 . The stoichiometry of passive transporters defines the directions that span , while active transporters define the displacement of 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 , the system evolves toward states of lower energy, but, as mentioned above, the state vector can only move along the two linearly independent directions and . In other words, is restricted to an affine plane containing . If we denote by × a normal vector to , the subspace of attainable states is characterized by the equation
where is the transpose of . 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 . Now, assume that the hydrogen pump ejects x moles of hydrogen before stopping, bringing the system to the state (see 2.18). To find an equilibrium, the state is now restricted to another subspace, denoted by , whose equation reads
It is notable that, for , the subspaces and are disjoint ( takes different values), and as a consequence, the resulting system equilibria are distinct.
Fig 4 shows a slice view of the space of states in the plane, with representations of the subspaces 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 is the point where a level curve is tangent to . As moves, the equilibrium state is pushed toward higher levels of energy, following the red curved trajectory.
Slice view of the space of states in the plane (n𝐇+, nCl−) for the toy problem.The level curves in the background represent the energy function G. Once the hydrogen pump has brought the system from 𝐧0 to 𝐧0+x𝐚, the passive evolution toward an equilibrium is only possible along the affine subspace 𝒱x, which is directed by 𝐯1 and 𝐯2 (𝐯2 is out of the plane). The minimum of G in 𝒱x corresponds to points where the level curves of G are tangent to 𝒱x. As x increases, the plane 𝒱x continuously moves, dragging the equilibrium state 𝐧(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.*
In the toy model from Fig 3, the extent of active transport x is the input variable, while the output is the corresponding equilibrium state . When the amount of hydrogen ejected by the pump describes an interval [0,x], it progressively displaces the passively attainable subspace , and the system equilibrium is continuously brought from to . A more general multi-membrane complex might feature several active transporters, whose chemical extents are stored in a vector , and whose corresponding directions are stored as the columns of a matrix . In that case, the system equilibrium is a function of the current extent of all active transporters, defined by
The space can always be characterized by an equation of the form
where the matrix has N rows and columns, and the columns of form a basis of .
A quasi-static simulation scenario. In this paper, the objective of the simulation is to evaluate the function when 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 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, simulation results look like an evolution as a function of a scalar variable. For instance, the example from Fig 1 involves two active transporters. If we run the simulation for with s describing the interval , 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 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, provided in S1 File, 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 evaluating 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 , the corresponding equilibrium is defined by
with the notation from (23) and (24). It can be shown that the solution to (25) satisfies some necessary first-order optimality conditions, also known as Karush-Kuhn-Tucker conditions [38, Section 12.3]. These conditions read
where is a Lagrange multiplier of size . In (26), the first line describes the balance of forces that characterizes the equilibrium point. The multiplier is adjusted so that the artificial force enforces the constraint . The second line corresponds to the constraint itself in matrix form.
To obtain a simpler system, we seek to eliminate . The matrix , whose columns form a basis of , does not depend on , as the are parallel to each other. We define the matrix , with N rows and columns, whose columns form a basis of . In particular, the columns of the block matrix form a basis of , and . By applying to the first line of (26), we obtain the equivalent system of size N
This system illustrates well that the coordinates of along directions collinear to are defined by the balance of forces, while the coordinates along directions orthogonal to are explicitly determined from .
Let us write (27) in the case of the toy problem from Fig 3, with and defined by
As the derivatives of and do not raise numerical issues, we only specify the expression for the chemical energy derivatives using (4) and (5), leading to the system
In practice, we aim to solve (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 (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 instead of functions of . Then, N new lines are added to (29) to enforce a relation of the form , where f and g are functions defined on and chosen to exhibit good numerical properties, namely
so that the enforced relation is either
depending on the values of and (note that may take values out of [0,1] if has negative components). As a consequence, the relaxed system, which is actually solved in the numerical procedure, reads
By solving (32) successively for x describing an interval 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 Fig 1 which symbolizes a guard cell.
3. Results
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 describe the parameters that we use to approximate the configuration of a guard cell, along with the choice of transporters (Sect 3.1). In Sects 3.2 and 3.3, we run simulations to explore the possibilities offered by the model. In Sects 3.4 and 3.5, we focus on a sensitivity analysis based on the second derivatives of the energy function. This analysis leads to a physical interpretation of the simulation results (Sect 3.6).
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 ( ), chloride ( ), potassium ( ), water ( ).
The guard cell is represented by the two-membrane complex summarized in Fig 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
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–46]. In the plasma membrane, we introduce a chloride/hydrogen symporter (which transports both species in the same direction), that has not yet been molecularly identified [40,47]. 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: + . Using different stoichiometric coefficients for the plasma membrane hydrogen/chloride symporter, such as or , 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.
Table 1: Numerical values of chosen parameters and initial conditions for the simulation.
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 for the Young modulus in the longitudinal direction and − for the wall thickness [22–24], which should result in a surface stiffness modulus × − . In the computational model, we set the surface stiffness modulus to , 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 Fig 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 , with x1 denoting the progress (in mol) of the plasma membrane pump, and x2 denoting the progress of the vacuole membrane pump. As mentioned in Sect 2.2, x1 and x2 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 between the two pump rates, i.e. . For each new value of s, the simulation computes a corresponding equilibrium state . Fig 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 across membranes, while the hydrostatic pressure and electric potentials across membranes are defined by
respectively (see (12) and (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 (given a molar volume of water of ), 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 − to . 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 ⋅ , while the amount of hydrogen actively ejected by the pump is 2 ⋅ .
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.
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 Fig 5, the ratio between pump rates is , 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 α, namely 40%, 80% and 100%. We do not explore values of α greater than 100%, as they raise numerical instabilities due to chemical amounts in the cytoplasm becoming close to zero. Results, gathered in Fig 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 Figs 5–7.
Simulation results with constant ratios between pump rates.The rates of the vacuole pump are α=40% (left), α=80% (middle), and α=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.
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.
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. x1(s) = s), and then we adjust the vacuole pump progress x2(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 Fig 7. The vacuole pump rate can be read as the slope of the vacuole pump progress in Fig 7b. For each value of the progress s, it is equal to the current fraction of cell volume occupied by the vacuole (Fig 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 Fig 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 (Fig 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. (Fig 5a). The massive intake of chloride is accompanied by an approximately equal intake of potassium (Fig 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 (Fig 5h). Finally, water intake (Fig 5b) is synonymous with volume increase, and the reaction of the elastic cell wall to this deformation can be seen in Fig 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 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 Sect 3.6.
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 [6,51,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 [6,51,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
where is a species able to bind with hydrogen. When the hydrogen concentration decreases, the reaction provides some hydrogen by consuming , while the reaction produces 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 in the cytoplasm. In the same fashion as passive transporters, this chemical reaction is represented by one of the vector directions that span , as it gives the system the freedom to evolve toward production or consumption of .
Fig 8 shows the simulation results for the buffered model, with the same pump rates as in Fig 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 Fig 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 Fig 8 show that our model can be extended to produce more realistic simulations.
Simulation result for the guard cell model with hydrogen buffer.The plotted physical quantities are the same as in Fig 5. The cytoplasm pH values follow a more realistic behavior when hydrogen buffering is taken into account.
3.4. Sensitivity analysis and second derivatives of the energy function
The next sections are dedicated to a sensitivity analysis of the guard cell model, based on second-order information about the energy function. Our aim here is to explain the ranges 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.
In this section, we emphasize through a short calculation the role of the energy second derivatives in determining the system response to a perturbation, and then we say a word about the physical meaning of the energy second derivatives. As G is a regular function of several variables , its second derivatives at are stored in a symmetric matrix called the Hessian matrix, denoted by , and the Hessian coefficient at index (i,j) is denoted by .
First, let us clearly identify where the energy second derivatives appear in the calculations. Remember that, for a fixed , the equilibrium state is solution to the nonlinear system (27). Now, assume that is subject to a perturbation . The equilibrium perturbation is the solution to a tangent linear system involving the Hessian matrix of the energy function. This linear system, obtained by differentiating (27), reads
Like previously, the first line of (36) concerns the coordinates of in the direction of , which adjust to re-establish a physical equilibrium, while the remaining coordinates are explicitly determined from by the second equation. The Hessian matrix plays a role in the evolution of in the directions of passive transporters, as a reaction to an imposed change in the directions orthogonal to . 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 extremity, causing it to extend by a length x in a one-dimensional way. The spring elastic energy function reads
Also, and 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 , the spring length changes by . 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 , and as a consequence, a force perturbation causes the length to change at first order by , where 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 . 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 , which makes it possible 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 xi, and due to the imposed total length, the spring deflections all add up to an imposed value y. The xi values are determined by solving the minimization problem
where ki is the stiffness of the i-th spring. Now, if the imposed cumulated deflection y changes by a perturbation , the system will find a new equilibrium state, where each deflection xi has changed by a small displacement
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 , then, solving a minimization problem similar to (38) yields the first-order responses to a perturbation ,
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 Sect 3.2.
3.5. 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.
Fig 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 , but their structures remain similar all along the simulation (see Hessian matrices evaluated at the initial state in S2 figS2 Fig). Concerning the axis labels, each row and column of a Hessian matrix is associated with one coordinate of . Chemical amounts are organized in the same order as in (33), 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 and . The values written in Fig 9 are in logarithmic scale, i.e. 27.2 on the matrix heatmap means that the corresponding coefficient is equal to 10^27.2^ in absolute value, while a blank box denotes a null coefficient. Therefore, a number visible in the total Hessian heatmap (Fig 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.
Hessian matrices for each term of the energy function, evaluated at the final state of the simulation.Hessian coefficients (in J/mol2) are represented in logarithmic scale, i.e. 27.2 means 1027.2 and a blank case means that the coefficient is zero. The Hessian matrices define the sensitivity of chemical, electrostatic and mechanical forces to perturbations.
To investigate the important balances in the system, we plotted in Fig 10a the eigenvectors , and eigenvalues of the total energy Hessian matrix, evaluated at the final state. Fig 10b shows the evolution of eigenvalues during the simulation. Except for , the eigenvalues remain well separated and keep the same ordering. In Fig 10c, we plot the angle between each eigenvector and its final direction . 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.
Spectral information about the total energy Hessian matrix.(a) Eigenvectors 𝐮1,…,𝐮8 of the total energy Hessian at the final state, with corresponding eigenvalues λ1,…,λ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 𝐮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.
Working with eigenvectors allows one to break the space of states into a set of orthogonal directions , each one representing a physical process. The system resists to displacements in each direction independently of the others, to an extent given by the eigenvalue . It is striking 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 and appear in the chemical energy Hessian matrix as diagonal parameters, as well as and . The eigenvectors and , correspond to variations of hydrogen amounts in the vacuole and cytoplasm, respectively, while and 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 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 and correspond to physical processes where chloride and potassium amounts evolve in opposite directions, resulting in charge transfers across membranes. Moving in the direction results in transferring charges from the cytoplasm to the vacuole and external environment, while in , charges are moved from both compartments to the external environment. Also, appears in the electrostatic Hessian matrix, and is close to certain coefficients in the same matrix (see framed zones in Fig 9b). The nondiagonal structure of the electrostatic Hessian matrix comes from the two terms of the electrostatic energy expression (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 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 in Fig 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 represents same-sign variations of water amounts in both compartments, resulting in variations of the cell volume. The corresponding eigenvalue 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 (S2 Fig a). The fact that 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 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 Fig 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 Fig 5. In Fig 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 ( ) are several orders of magnitude smaller than variations of chloride and potassium ( ), 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 Fig 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 ( ) in the system, we plotted in Fig 5h the charge imbalance across each membrane, defined by (see 2.8). It denotes the amount of charges (in mol) stored in each membrane/capacitor (see the circuit analogy in Fig 2). Like previously, the annotations show the charge imbalance variations at the plasma membrane ( ) and at the vacuole membrane ( ). Because charge imbalance shares units with chemical amounts, its variations (Fig 5h) can be compared with the amount variations of all species (Fig 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 ( ). 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 ( ). 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 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 Fig 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.
3.6. 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 through the prism of this hierarchy between the terms of the energy function, going from the activation of hydrogen pumps all the way to water entering the cell. Our simulation results in Sect 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 Fig 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 Fig 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.
Summary of the transport cascade at the plasma membrane in the guard cell model without buffering, showing the hierachy 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.
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.
4. Discussion
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 resulting 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 [6,9,28] 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 Fig 4, a necessary condition for to change as a function of x is that the subspaces and be distinct for two values , i.e. the passive transporters do not allow the system to come back to its initial state . In contrast, if the direction is collinear to in the figure, then is an element of for every x in . Therefore, the energy minimizer in always remains , 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 is not a linear combination of the passive directions .
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 dynamic model using the energetic approach is given in [9]. With our notation, the flux J across a passive transporter associated with the direction is proportional to the energy derivative along , i.e. , where α is a conductance parameter specific to the considered transporter. 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 conductance. However, it can be used to create more complex behaviors, for instance by considering transporters whose dynamics is affected by voltage or pH.
Supporting information
S1 FileCode used to run simulations and generate figures.The archive contains the Python code to reproduce simulations presented in this paper. Instructions to run the simulations are given in the file Readme.md.(ZIP)
S2 FigHessian matrices at the initial state of the simulation.This figure is similar to Fig 9, but the Hessian matrices are evaluated at the initial state of the simulation . In particular, the analysis of Hessian matrices from Fig 9 still apply to these matrices, as similar differences between coefficient orders of magnitude are present here.(EPS)
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Martinoia E, Meyer S, De Angeli A, Nagy R. Vacuolar transporters in their physiological context. Annu Rev Plant Biol. 2012;63:183–213. doi: 10.1146/annurev-arplant-042811-105608 22404463 · doi ↗ · pubmed ↗
- 2Krüger F, Schumacher K. Pumping up the volume - vacuole biogenesis in Arabidopsis thaliana. Semin Cell Dev Biol. 2018;80:106–12. doi: 10.1016/j.semcdb.2017.07.008 28694113 · doi ↗ · pubmed ↗
- 3Kaiser S, Scheuring D. To lead or to follow: contribution of the plant vacuole to cell growth. Front Plant Sci. 2020;11:553. doi: 10.3389/fpls.2020.00553 32457785 PMC 7227418 · doi ↗ · pubmed ↗
- 4Dünser K, Schöller M, Rößling A-K, Löfke C, Xiao N, PařízkováB, et al. Endocytic trafficking promotes vacuolar enlargements for fast cell expansion rates in plants. Elife. 2022;11:e 75945. doi: 10.7554/e Life.75945 35686734 PMC 9187339 · doi ↗ · pubmed ↗
- 5Vitali V, Sutka M, Amodeo G, Chara O, Ozu M. The water to solute permeability ratio governs the osmotic volume dynamics in beetroot vacuoles. Front Plant Sci. 2016;7:1388. doi: 10.3389/fpls.2016.01388 27695468 PMC 5024706 · doi ↗ · pubmed ↗
- 6Li K, Grauschopf C, Hedrich R, Dreyer I, Konrad KR. K+ and p H homeostasis in plant cells is controlled by a synchronized K+ /H+ antiport at the plasma and vacuolar membrane. New Phytol. 2024;241(4):1525–42. doi: 10.1111/nph.19436 38017688 · doi ↗ · pubmed ↗
- 7Kahm M, Navarrete C, Llopis-Torregrosa V, Herrera R, Barreto L, Yenush L, et al. Potassium starvation in yeast: mechanisms of homeostasis revealed by mathematical modeling. P Lo S Comput Biol. 2012;8(6):e 1002548. doi: 10.1371/journal.pcbi.1002548 22737060 PMC 3380843 · doi ↗ · pubmed ↗
- 8Astaburuaga R, Quintanar Haro OD, Stauber T, Relógio A. A Mathematical model of lysosomal ion homeostasis points to differential effects of Cl- transport in Ca 2+ dynamics. Cells. 2019;8(10):1263. doi: 10.3390/cells 8101263 31623161 PMC 6848924 · doi ↗ · pubmed ↗