## Abstract

An elastic membrane stretched between two walls takes a shape defined by its length and the volume of fluid it encloses. Many biological structures, such as cells, mitochondria and coiled DNA, have fine internal structure in which a membrane (or elastic member) is geometrically ‘confined’ by another object. Here, the two-dimensional shape of an elastic membrane in a ‘confining’ box is studied by introducing a repulsive confinement pressure that prevents the membrane from intersecting the wall. The stage is set by contrasting confined and unconfined solutions. Continuation methods are then used to compute response diagrams, from which we identify the particular membrane mechanics that generate mitochondria-like shapes. Large confinement pressures yield complex response diagrams with secondary bifurcations and multiple turning points where modal identities may change. Regions in parameter space where such behaviour occurs are then mapped.

## 1. Introduction

The mitochondrion is an organelle of the cell that is responsible for converting sugars into adenosine triphosphate (ATP) and is for good reason considered the ‘engine’ of the cell [1]. This double-membrane structure, shown in figure 1, consists of a permeable outer membrane, which allows transport of sugars into the intermembrane space, bounding an impermeable inner membrane where the production of ATP takes place. The efficiency with which a mitochondrion processes sugars is directly related to the surface area of the inner membrane, which can be anywhere from 4 to 10 times as large as the surface area of the outer membrane. The typical inner membrane has size 1 µm and is characterized by a bending rigidity *α* = 10^{−12} erg, tension *β* = 10^{−3} ∼ 10^{−1} pN nm^{−1} and internal pressure *p* = 4 × 10^{−2} atm [2,3]. To accommodate such large surface areas, the inner membrane deforms itself into a series of sharp folds (invaginations) known as cristae. In addition to energy production, the structure of the inner membrane can affect specific chemical pathways and processes, such as cellular differentiation and the cell cycle, through its geometry. The mechanics governing the morphology of the inner membrane is our focus.

As shown in figure 1, the inner membrane behaves like a closed elastic sheet that deforms as a two-dimensional object. This sheet fits within the outer membrane and has three-dimensional aspects, but it is expected that the mechanics should be described well with a two-dimensional model. In this study, we develop a two-dimensional Cartesian model of a mitochondrion where we hold the outer membrane as a rigid confining surface and focus on the shape of the inner membrane. The inner membrane is characterized by a tension *β* that encloses a volume of fluid *V* (intermembrane space) with pressure *p*. The interactions between the inner membrane and outer membrane are modelled using a general confining pressure that prevents the inner membrane from penetrating the outer surface. Our idealized model for the observed confinement (outer membrane) geometry should reproduce the essential feature of the interaction, especially in the limiting case having a large number of folds.

Biological interfaces/membranes are subject to instabilities such as surface wrinkling, folding and creasing in the presence of (i) geometrical constraints (e.g. particles [4,5], fibres [6,7] and thin films [8]) or (ii) external loadings (e.g. mechanical forces [9], hydration stresses [10–12] or capillarity [13,14]). The review by Li *et al*. [15] catalogues the complex morphology of observed ‘surface instabilities’ that occur in soft biological materials and illustrates the relevance of such phenomena to industrial applications, such as soft lithography, flexible electronics and biomedical engineering. The transition from wrinkle to fold is seen upon increasing the membrane compression to a given threshold [16]. In unconfined geometries, the shape of surface wrinkles may be computed using any number of mathematical techniques, such as elastic rod theory or by considering a fluid–structure interaction problem of a periodic, length-preserving bilipid membrane modelled by the Helfrich energy immersed in a viscous fluid [17]. The latter technique has been used to simulate the flow of red blood cells [18]. Here, we use numerical continuation to compute membrane shapes [19].

Continuation methods allow one to ‘follow’ a solution to a given problem as a parameter changes. The method is well known in a number of different fields: astrophysics [20–22], vortex flows [23], capillary surfaces [24] and rigid body dynamics [25], but has yet to become popularized with elastic membranes to our knowledge. We use the method to mimic the mechanics of the inner membrane of the mitochondrion. In essence, we perform ‘numerical experiments’ and compare our predictions against experimental observations. We consider two scenarios. The first treats the mitochondrion as a ‘closed system’ with fixed fluid volume. We then calculate the membrane response (*p*, *β*) to increased length *S*. The second assumes the membrane tension is a material parameter (fixed) and the system is ‘open’ to allow fluid transport through the ‘intermembrane space’. It is clear from the comparison of our results that the second scenario most probably governs the mechanics of the inner membrane, whereas the first scenario is potentially relevant to other applications such as wrinkling [26] and folding [27,28] in elastic sheets.

From a broader mechanical perspective, the shape of elastic membranes belongs to a class of problems classically referred to as an ‘elastica’. Examples include, but are not limited to, vesicles whose shapes can be determined using analytical techniques [29,30] or numerical techniques such as the finite-element methods [31] or the discrete space variation method [32]. More recently, attention has been paid to solid elastic members geometrically constrained to lie within a rectangular [33] or circular channel [34]. Typical modelling of such situations involves introducing geometric pseudo-forces to account for the presence of the geometrical support [9,35]. Alternatively, one can introduce intermolecular forces, as done with the adhesion of elastic sheets [36], snap-through of graphene sheets onto complex substrate geometries [37], peeling of elastic sheets [38] and the pull-in instability of carbon nanotube tweezers [39].

An elastic member in contact with a fluid reservoir can deform from the associated fluid pressure, either hydrostatic or capillary. For hydrostatic pressures, the buckled shapes of thin elastic sheets have a closed-form solution [40,41] with post-buckling solutions numerically computed using continuation techniques [42]. For liquid interfaces endowed with surface tension, e.g. drops, capillary pressure and other wetting forces [43] give rise to ‘elastocapillary deformations’ with applications such as capillary origami [44], wrapping of droplets by elastic sheets [45] and the formation of wetting ridges [46].

The construction of the governing equations through a variational principle subject to auxiliary constraints has been well studied in related problems on capillary surfaces (membranes with negligible elastic resistance), whereby the stability of solutions to the Euler–Lagrange equations associated with the respective energy functional can be read off a preferred bifurcation diagram [24,47]. Although we do not discuss stability of solutions here, the relationship between these mathematical theorems and stability of the shapes predicted in this manuscript would certainly be of interest.

We begin this paper by formulating a mathematical model that describes the shape of an elastic membrane confined by two parallel surfaces by introducing a confinement pressure. Families of membrane shapes are computed according to two alternative representations for the mechanics of the inner membrane of the mitochondria. Response diagrams which plot the membrane tension and fluid pressure as functions of total arclength are computed. For large confinement pressures, the membrane response often exhibits secondary bifurcations (SBs) and multiple turning points where the membrane's modal identity may change. We catalogue the canonical behaviours and map the corresponding regions of parameter space. Our results predict membrane shapes close to observed morphologies of the inner membrane and provide insight into the governing mechanics. Finally, some concluding remarks are offered regarding the relevance of our model to experimental observation and the applicability of our methods to related problems.

## 2. Mathematical formulation

Consider an elastic membrane ∂*D* of length *S* held by tension *β* that encloses a domain *D* with volume of fluid per unit length *V* and pressure *p*, as shown in figure 2. The membrane shape *h*(*x*) is defined in Cartesian coordinates and ‘confined’ to lie within the domain bound by two rigid walls at *y* = ±*H*, with incompressible fluid of fixed volume above and below the interface. Our model can be viewed as an idealized mitochondrion by treating the solid wall as the outer membrane, the volume of fluid as the intermembrane space and the interface *h* as the inner membrane. The energetics of this configuration are described by the Helfrich [48,49] free energy functional
2.1where *κ* and *α* are the curvature and bending rigidity of the membrane, respectively, and *p*_{c} is a confinement pressure introduced to model the interactions between the membrane and bounding surfaces. Although referred to as tension herein, note that *β* > 0 corresponds to compressive forces. We are interested in membrane shapes on a periodic domain and therefore introduce a horizontal length scale *L*. Lengths are then scaled by *L*, while tension *β* and pressure *p* are scaled by the bending rigidity *α*;
2.2

Herein, we drop the hats for simplicity.

The Euler–Lagrange equations for the functional (2.1) are given by
2.3where Δ is the surface Laplacian [50] and we have introduced a general confinement pressure *p*_{c}(*h*). We model the interactions between the membrane and the bounding solid using the power-law form,
2.4whose magnitude is controlled by the parameter *M*. Here 2*H* is the wall spacing that defines confinement geometry, whereas the exponent *n* models the interaction, e.g. *n* = 3 for medium range van der Waals [36,37,51]. For simplicity, we set *n* = 1 and *H* = 1/2 unless otherwise stated. We choose the power-law form (2.4), because it is a ‘hard’ potential that prevents the membrane from penetrating the solid. It is straightforward to analyse other functional forms for the repulsive confinement pressure, such as the ‘soft’ potential for hydration stresses [12,52] (see appendix A).

The membrane equation (2.3) is coupled to an equation for the curvature *κ* = *κ*(*h*), written in Cartesian coordinates, which will be defined shortly. These governing equations are augmented by the following boundary conditions:
2.5which are periodic conditions (with period 2*L*), where we have taken the symmetric extension about the midplane *x* = 0. Note that this is tantamount to setting the phase of the membrane shape. Finally, the integral conditions
2.6

fix the volume *V* of fluid and total arclength *S* of the membrane, respectively. Note that these are defined on the half-domain.

## 3. Results

The energy functional (2.1) can be viewed as a constrained variational principle provided one treats the tension *β* and pressure *p* as Lagrange multipliers (unknowns) that are fixed by the integral constraints on volume *V* and arclength *S*. Here, (*β, S*) and (*p*, *A*) are dual variables. One can consider the tension *β* and pressure *p* as known parameters, in which case *V* and *S* are solution measures. This dual construction lends itself to computational methods that involve continuation, whereby a known solution is ‘continued’ in one or more parameters to obtain families of solutions [53]. We use the continuation software package AUTO [19] to compute families of membrane shapes. With regard to the morphology of the inner membrane of the mitochondria, the strategy we employ is to perform ‘numerical experiments’ using continuation that represents the mechanics typical to the inner membrane.

### 3.1. Long-wave approximation

We begin by examining the long-wave limit with curvature given by . In this limit, equation (2.3) becomes 3.1

For fixed *V*, the ‘unconfined’ (*M* = 0) solutions,
3.2belong to the one-parameter family shown in the bifurcation diagram of figure 3*a*, which plots tension *β* against arclength *S*. Solution branches emanating from the trivial branch *S* = 1 are characterized by a wavenumber *k* and parametrized by arclength, as shown for the *k* = 1 solution in figure 3*b*. Note that the membrane tension is constant, regardless of the arclength. The response diagram illustrates the dual construction mentioned above: (i) traversing the graph horizontally for prescribed arclength *S* reveals the computed tension *β* (Lagrange multiplier), while (ii) traversing the graph vertically for fixed tension gives the arclength of that solution.

The response of a membrane in confinement (*M* ≠ 0) becomes more complicated, as shown in figure 4*a* for the branch *k* = 1. Here, we compute the response diagram for fixed *V* = 0.5 by continuing in *S* from the trivial state *h* = 0 and contrasting against the unconfined (*M* = 0) response. Both confined and unconfined branches bifurcate from the trivial branch at *β* = *π*^{2} and, upon increasing *S*, the membrane in confinement responds by increasing its tension (stiffening) in order to accommodate the additional membrane length. The corresponding membrane shape tends to steepen on the sides and hug the bounding surface when compared with the limiting unconfined solution that obeys the geometrical constraint, as shown in figure 4*b*.

### 3.2. Full model

The governing membrane equation (2.3) is coupled to the nonlinear curvature 3.3

Recall that we have used Cartesian coordinates for convenience to simply analyse the confinement pressure and avoid self-intersecting shapes. This construction ensures single-valued functions *h*(*x*). As mentioned above, we use continuation methods to solve the nonlinear problem in order to better understand the mechanics of the inner membrane by comparing our predictions against observed morphologies. We explore two possible scenarios: (i) fixed *V* and (ii) fixed *β*. In the first case, we fix *V* and use *S* as the principal continuation parameter, adjusting the tension and pressure accordingly. That is, we study how the membrane responds (*p*, *β*) as more arclength is packed between the two ‘outer membranes’. By contrast, for fixed *β* we vary *p* and observe *V* and *S* for each (*p*, *β*). The two alternatives highlight the dual construction of the variational principle.

#### 3.2.1. Prescribed fluid volume

Figure 5*a* plots the response diagram *β*−*S* for fixed *V* = 0.5 with confining pressure magnitude *M* = 0.1. Note that the membrane initially responds by slightly decreasing its tension (softening) as arclength increases, which is different from the long-wave model of figure 4*a*. We attribute this feature to the nonlinear terms in the model which allows for high local stress to be balanced by the confinement pressure terms instead of surface tension. As arclength increases further the membrane begins to feel the bounding surface and responds by increasing its tension. For the mode *k* = 2, the tension dramatically increases from 25 to 160 despite a marginal increase in arclength. Limiting cases of membrane shapes with the largest arclength on the respective branches are shown in figure 5*b*. Our computational approach did not find solutions past these values of *β*, which is not surprising given the large curvatures and steep sides observed in figure 5*b*. Note the steep sides reminiscent of the inner membrane shapes in figure 1.

Continuation methods also allow one to ‘continue’ different solution measures, such as the bifurcation points shown in figure 5*a*, in a given parameter. Figure 6 plots the locus of bifurcation points in different parameters; (*a*) confinement pressure *M*, (*b*) *V* and (*c*) *p*. Interestingly, for large confinement pressure *M* the order in which membrane shapes appear can become distorted. For example, traversing the graph by increasing tension *β* for fixed *M* = 1 intersects (in order) the *k* = 1, 2, 3 branches, whereas for *M* = 80 the order becomes distorted *k* = 2, 1, 3; the *k* = 2 shape appears before the *k* = 1 shape. We attribute this feature to the strong interactions with the solid surface through the nonlinearity of the confinement pressure.

Along the same lines, the locus of bifurcation points in the volume *V* space reveals that starting from the neutral state *V* = 0.5 the membrane dramatically increases its tension and pressure as *V* is removed. Note the symmetry for increasing volume *V* from the neutral state. With regard to the mitochondrion, large pressures and tensions may be expected to result in rupture of the inner membrane under such circumstances. As such, figure 6*b* could be viewed as a guide in controlling the hydration (equivalently, *V*) or *p* to ensure the membrane tension remains below a critical rupture threshold. As with the confinement pressure *M*, spectral reordering may occur when the membrane strongly feels the presence of the confining surface.

In general, the membrane response becomes more complex with increasing *M*, as shown in figure 7*a* which plots *β* against *S* for *M* = 1. Initially, for increasing *S*, the membrane responds by softening (*β* decreasing) until the first limit point (LP), defined by the appearance of two solutions, is traversed, after which the membrane stiffens (*β* increasing) until the second LP, whereby the membrane changes its modal identity from *k* = 1 to 3, after which it continues to soften. Our continuation method allows us to follow the second LP, thereby constructing the curve in the *β*−*M* parameter space where the modal identity changes. The locus of LPs shown in figure 7*b* is a non-monotonic curve that could not have been predicted either *a priori* or by long-wave theory. Needless to say, cataloging membrane shapes for large *M* becomes increasingly difficult.

SBs off the primary branch are typical for larger values of *M*. Figure 8*a* plots the response diagram that exhibits an SB off the primary branch *k* = 2 for *M* = 0.25 and *V* = 0.5. Solutions on the secondary branch preserve the modal identity *k* = 2, but exhibit a different symmetry. Both solutions are symmetric about the centre line *x* = 0. However, solutions on the primary branch are also symmetric about the vertical generator *x* = 1/2, while solutions on the secondary branch are not. The locus of SB points for the branch *k* = 2 are shown in figure 8*b* in the *β*−*M* parameter space. As could be expected, larger confinement pressures *M* lead to increased membrane tension *β*. The classification of the space of membrane shapes can become quite complex, as there are three equivalent membrane shapes for *β* = 30. The question of which membrane shape occurs in the mitochondrion may depend upon the loading path and can only be answered by stability analysis, which is outside the scope of this paper.

#### 3.2.2. Prescribed membrane tension

The second scenario we consider is one where the membrane tension *β* is fixed and the pressure *p* is allowed to vary. As stated earlier, in this case the volume and arclength are solution measures, *β, p* are known material parameters. The response diagram of figure 9*a* plots the pressure *p* against arclength *S* for fixed *β* = 10. As shown, the pressure dramatically increases with the arclength due to the strong interactions between the membrane and confining surface; a large fraction of the membrane shape lies near the idealized ‘outer membrane’. The limiting case corresponds to a total arclength *S* ≈ 2.84, which is nearly 95% of the length *S* = 3 of the confining slot geometry with *H* = 0.5. Figure 9*b* plots the volume *V* against arclength *S* to show volume is removed as pressure increases. The limiting case has volume *V* ≈ 0.058. The morphology of the membrane shape shown in figure 9*a* (inset) closely resembles the inner membrane of the mitochondrion observed in the literature (figure 1), leading us to believe that one should treat the membrane tension as a *known* material parameter. In addition, our computed shapes qualitatively resemble those computed by Guo & Li [32], who introduce a repulsive intermembrane potential (disjoining pressure) to model self-avoiding behaviour and prevent membrane self-intersection. Note that our shapes are non-intersecting via construction in Cartesian coordinates.

### 3.3. Biological interpretation

We revisit our results to interpret the biological mechanics of the inner membrane of the mitochondrion. The first scenario could be idealized as a ‘closed system’ in which the fluid volume *V* is fixed and the membrane responds (*p*, *β*) to increasing length *S*. Note that our computed membrane tensions *β* in figures 5*a*, 6 and 8 are all consistent with observed values for the mitochondrion (see the Introduction) [2,3]. Figure 7 shows the response can be complex; the membrane initially softens, followed by a period of stiffening and then softens again after changing its wavenumber character. In this case, the membrane tension must actively change throughout the process and this must be accomplished through the chemistry of the lipid bilayer. Although it is well known that biochemical processes can be complex, this scenario seems unlikely.

By contrast, the second scenario, where the membrane tension is treated as a material parameter, can be idealized as an ‘open system’ that allows for fluid transport through the ‘intermembrane space’. Our results predict membrane shapes that closely resemble the morphology of the mitochondrion (figure 9). Additionally, the limiting case has a large pressure *p* ∼ 10^{4} whose order is consistent with experimental observations [2,3] ( calculated from the numerical values given in the Introduction). These two observations strongly suggest that the inner membrane obeys the mechanics of the second scenario.

## 4. Concluding remarks

We have developed an idealized model of an elastic membrane in confinement in order gain insight into the mechanics of the inner membrane of the mitochondrion, a double-membrane structure, whose outer membrane acts as a confining surface. The governing equations are derived from a constrained variational principle that has a dual interpretation that we exploit to perform ‘numerical experiments’ that replicate typical mechanics that may lead to the unique morphology of the inner membrane. Recall that we have chosen to use a Cartesian representation to avoid self-intersecting shapes, which restricts the class of allowable functions. For example, shapes with a bulbous tip are not allowed here but would be in an arclength parametrization. This construction could affect the response diagrams (e.g. figures 5*a* and 7*a*). We compute the bifurcation diagrams using the continuation software AUTO to show how the membrane's tension and associated fluid pressure *p* respond to increases in total arclength. For small confinement pressures, the membrane initially softens as its length grows until it feels the presence of the confining surface, after which the membrane stiffens significantly in order to accommodate additional length (figure 5*a*). As the confinement pressure *M* increases, the membrane response can become complex, often exhibiting multiple turning points and SBs (cf. figures 7 and 8). Our work has catalogued such behaviour and mapped corresponding regions in parameter space.

Our results indicate that treating the membrane tension *β* as a material parameter in an ‘open system’ that allows for fluid transport through the ‘intermembrane space’ yields membrane shapes that compare favourably with observed morphologies, thereby generating insight into the mechanics of the inner membrane. This interpretation is consistent with the mechanism by which the mitochondrion controls the transport of sugars across the outer membrane into the intermembrane space to be processed by the inner membrane. Clearly, it is crucial to include the fluid mechanics of the transport process in developing future models of the mitochondrion. For example, how does dehydration or rehydration affect the efficiency with which the mitochondrion can process sugars into ATP?

From a mathematical perspective, we hope our results show that the application of continuation methods, which are typically applied to problems on capillary surfaces, can be used on related problems involving elastic membranes. Although we do not focus on stability issues here, the constrained variational structures inherent in problems on elastic membranes are well suited to bifurcation theoretic approaches, such as the Poincaré–Maddocks turning point theorems [54].

Extensions to our model will involve changing the geometry of the confining surface to better reflect the outer membrane of the mitochondria or increasing the degree of freedom in the model by treating the confining surface as a deformable membrane, whose shape must also be determined as part of the problem solution. Along a similar thread, there are many problems in biological mechanics that involve biological members ‘confined’ by some geometry, such as DNA coiling, that may be applicable to our general solution procedure. Perhaps the most direct extension would be to consider the three-dimensional case; this is non-trivial and involves solving a complicated partial differential equation associated with the free boundary value problem.

## Authors' contributions

J.B.B., M.J.M. and S.H.D. conceived the study and wrote and revised the manuscript. J.B.B. performed the analysis and simulations. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

M.J.M. acknowledges partial support from NSF grant no. DMS-1312935.

## Acknowledgements

We wish to thank Alejandro Sanchez, who did some of the initial calculations presented in appendix A. We also would like to thank Jeff McFadden for several helpful conversations.

## Appendix A. Hydration stresses

It is straightforward to analyse a membrane in the presence of a confinement pressure resulting from hydration stresses [12,52],
A 1by substituting into equation (2.3). Here *λ* is a characteristic length scale over which the hydration forces act. Repulsive hydration forces have been observed experimentally by Israelachvili & McGuiggan [55] and Leikin *et al*. [56]. Equation (A 1) can be viewed as a ‘soft’ potential whose magnitude remains finite as the membrane touches the solid support. By contrast, the hard potential (2.4) used in the main text gives rise to infinite forces as the membrane approaches the solid that prevents penetration into the solid.

Figure 10 plots the bifurcation diagram *S* against *β* in the long-wave limit for a membrane in the presence of hydration stresses. Note that the response diagram is qualitatively similar to that reported in figure 7 for the hard potential (2.4), with the exception that the membrane intersects the solid for these values of the parameters *M*, *λ* (see inset figure 10). That is, these computed membrane shapes do not physically represent the observed mitochondrion shapes that we seek to model. For the soft potential (A 1), there is a pseudo-support at some distance away from the actual solid, whose position cannot be predicted *a priori*. In figure 11, we plot the membrane shape for fixed arclength *S* = 3, as a function of *λ*, to show that the location of the pseudo-support depends strongly on (i) *λ* and (ii) the wavenumber *k*. For fixed *λ*, the *k* = 2 shapes may not intersect the solid while the *k* = 1 shapes may do so. Whether or not this occurs is not predictable *a priori* and it is primarily for this reason that we consider the hard potential (2.4) in this paper.

- Received May 25, 2016.
- Accepted June 24, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.