Compression, crumpling and collapse of spherical shells and capsules

The deformation of thin spherical shells by applying an external pressure or by reducing the volume is studied by computer simulations and scaling arguments. The shape of the deformed shells depends on the deformation rate, the reduced volume V/V0 and the Föppl–von Kármán number γ. For slow deformations the shell attains its ground state, a shell with a single indentation, whereas for large deformation rates the shell appears crumpled with many indentations. The rim of the single indentation undergoes a shape transition from smooth to polygonal for γ≃7000(ΔV/V0)− 3/4. For the smooth rim the elastic energy scales like γ1/4 whereas for the polygonal indentation we find a much smaller exponent, even smaller than the exponent 1/6 that is predicted for stretching ridges. The relaxation of a shell with multiple indentations towards the ground state follows an Ostwald ripening type of pathway and depends on the compression rate and on the Föppl–von Kármán number. The number of indentations decreases as a power law with time t following Nind∼t− 0.375 for γ=8×103 and γ=8×104 whereas for γ=8×105 the relaxation time is longer than the simulation time.


Introduction
A cell is the structural and functional unit of living organisms with a typical size of about 10 µm. There are many different cell types depending on the type of tissue that they are part of. Nevertheless, what all cells have in common is that they are the envelop of a portion of a complex fluid, the cytoplasm, containing many different kinds of small molecules, macromolecules and organelles in a highly organized and structured fashion. The envelop or cell membrane is in general a lipid bilayer with other molecules such as cholesterol and different proteins embedded in it. The cell membrane is fluid-like, that is, molecules can diffuse on the surface. Often the cell membrane is reinforced with a structure reducing the fluidity and introducing a shear modulus. For red blood cells this structure is the spectrin network [1]- [3]. Moreover, most cells have a dynamic cytoskeleton that maintains the cell's shape and plays an important role in cellular motion and transport within the cell [4]. In the process of endoand exocytosis, clathrin molecules assemble on the cell membrane to form a hexagonal pattern that induces the formation of a coated pit followed by a bud [5,6]. Apart from biological cells, spherical (icosahedral) viruses are another example of biological nano-cages with a shear modulus [7]- [10].
In recent years there have been serious efforts to synthesize artificial cells consisting of a lipid bilayer and an actin cytoskeleton at the micrometre level and to study their adsorption on surfaces as well as their mechanical properties [11]- [18]. Polyelectrolyte capsules with radii between 2 and 4.2 µm and a wall thickness of around 10 nm were synthesized by Gao et al [12]. Furthermore, the response of the capsules was measured after placing them in polymer solutions of various concentrations. In the osmotic deswelling experiments, the formation of single indentations was observed, and by fitting the critical osmotic pressure to the theoretical prediction from continuum elasticity theory an elastic modulus of about 500 MPa was extracted for these shells [12]. This is a value also found for plastics such as teflon or polypropylene. Similar-sized polyelectrolyte capsules were synthesized by Vinogradova et al, but in their experiments much more irregular shapes of osmotically deswollen capsules were observed [16]. Estimates of the Young modulus from antiferromagnetic (AFM) measurements were in the range of 70-250 MPa. For the osmotic deswelling of so-called spiderbags, elastic vesicles in the size range of 1-30 µm made of polymerized spidersilk electronmicrographs of the collapsed state (with a diameter of about 18 µm) appear highly crumpled [17]. These experimental results demonstrate that there is quite a variation in the shapes of collapsed shells. Similar crumpled shapes were found by Ratanabanangkoon et al [19,20] when giant lipid bilayer vesicles coated with streptavidin, which gives the vesicles a shear modulus, were slowly deflated using a micropipette. In applications, the permeability of the capsule wall is very important, because it is essential for keeping cargo inside or releasing it on demand. For this reason, capsules with different permeabilities have been constructed: some have very low permeability [19,21], while others have controllable permeability [22]. Therefore, it is important to understand capsule behaviour at either constant volume or constant osmotic pressure difference.
The systems discussed above have important biological functions or potential medical and technological applications. A detailed understanding of their mechanical properties requires a careful analysis based on statistical physics. Moreover, it turns out that there are great similarities between the structures formed in microscopic systems and in macroscopic systems such as an uninflated soccer ball or a table tennis ball. For these macroscopic systems we know that a soccer ball responds differently from a table tennis ball when pressing on it. The soccer ball is relatively soft and the rim of the indentation is smooth, whereas the table tennis ball is hard and forms an irregular indentation with the shape of a polygon. The same shape as the compressed table tennis ball was also observed in dried hollow silicon microcapsules [23,24].
For understanding the deformation of elastic membranes and shells, the Föppl-von Kármán number plays an important role, where R 0 is the radius of the shell and ν the Poisson ratio. For homogeneous elastic shells of thickness h, the two-dimensional (2D) Young modulus is given by K 0 = Y 0 h, with Y 0 the 3D Young modulus, and the bending rigidity by κ = Y 0 h 3 . An increase of γ leads to quantitative and qualitative structural changes of compressed shells. In particular, the formation of folds or stretching ridges, which are also characteristic for crumpled sheets, is a direct consequence of the inextensibility due to a large stretching modulus [25,26] and originates from the balance between stretching and bending energy. This can be readily seen when the progressive deformation of a flat thin elastic sheet is investigated. When such a sheet is deformed, for instance by confinement in a sphere, the first mode of deformation is the formation of a so-called d(evelopable)-cone that has zero Gaussian curvature everywhere [26]- [29]. Such a shape is formed by bending the elastic sheet without stretching it. However, when the confinement is too strong, the stress condenses in a disordered pattern of nuclei, each of them growing into a cone-like structure [29]. These structures interact with each other by forming stretching ridges or folds between each other [26]. The formation of the stretching ridges imposes a nonzero Gaussian curvature and requires the surface to be locally stretched in addition to the bending deformation. For homogeneous spherical elastic shells, the situation is fundamentally different and any type of deformation requires bending and stretching of the surface due to the coupling of bending and stretching in spherical topologies. The pressureinduced collapse of thin spherical shells is a very old problem in elasticity theory, and it has been known for a long time that for pressures larger than the critical pressure the spherical shape is unstable and collapses into a new shape with one or more indentations [30]- [32]. For subcritical pressures, the shell remains spherical but the surface is compressed. For such a deformation the bending energy is constant and the stretching energy rises rapidly. At the transition point it becomes more favourable to make indentations. This increases the bending energy but releases stretching energy. In figure 1 the characteristic shape and length scales involved in a single indentation are illustrated. The indentations are characterized by a depth H whereas the width of the rim d, is independent of the depth of the indentation. The radius of the indentation is geometrically determined and only depends on the size of the shell and the depth of the indentation (or the enclosed volume). This means that the elasticity of the surface imposes only a single length scale on the deformation of the surface. In this paper we investigate the different shapes that appear when a spherical shell is collapsed either under an external pressure or by a volume constraint. To illustrate this we present in figure 2 two representative examples of configurations that we find in the simulations. The reduced volume of the shells is in both cases the same but the Föppl-von Kármán numbers are different. The localization of the elastic energy in the highly deformed regions is displayed by the darker shades of grey. The shell with small γ (figure 2(A)) is characterized by a single indentation and looks rather regular whereas the shell with large γ (figure 2(B)) has several indentations and appears crumpled.
A simple scaling argument based on the balance of bending and stretching energy shows that for a single very small deformation the energy scales quadratically with the depth H of the indentation and crosses over to an H 3/2 scaling for larger indentations [32]. As long as indentations are non-interacting, i.e. when they are far apart on the surface, the total energy is just the sum of the energies of the individual indentations. We will show in section 3.2 that the ground state is always a configuration with a single indentation. Consequently, whenever configurations with multiple indentations are found the system is not in equilibrium. These nonequilibrium structures are mainly found, for large Föppl-von Kármán numbers, when stretching ridges are present. This paper is concerned with a very simplified model for elastic vesicles, capsules or cells, namely an elastic shell. An elastic shell is most easily thought of as a polymerized membrane. Due to their internal connectivity, these membranes have a shear modulus. The remainder of the paper is divided into two sections. In section 2 we briefly describe our simulation model and in section 3 we present and discuss the results of our simulations.

Simulation model
Our simulation model is based on a triangulated surface model for membranes and capsules [33]- [35]. A capsule of radius R 0 is modelled as a triangulated surface consisting of N v vertices located on the surface of a sphere. With a regular, icosahedral triangulation of the sphere, it is difficult to avoid the fact that the membrane parts near the 12 fivefold coordinated vertices with smaller triangles-corresponding to the corners of the icosahedronhave somewhat different elastic properties than those corresponding to the faces of the icosahedron (with larger triangles). Therefore, we employ a random triangulation. Such a triangulation is created by randomly distributing vertices on the surface of a sphere and equilibrating this configuration in a Metropolis Monte Carlo (MC) simulation. During the equilibration, vertices interact with the repulsive part of a Lennard-Jones potential. After equilibration, the configuration is 'frozen' and the neighbours of all vertices are determined using a Delauny triangulation [36]. This configuration defines the connectivity of the network. Clearly, one can obtain other different microscopic realizations of the system by doing more MC iterations. The spherical start configuration and neighbour list are used in further MC simulations.
In subsequent simulations, the neighbouring vertices are permanently linked via a harmonic potential that gives rise to a total stretching energy where the sum runs over all distinct nearest neighbour pairs of vertices i and j, r i j is the distance between two neighbours and r 0 i j is the equilibrium separation of a spring between vertices i and j. The rest length r 0 i j is determined at the start of the simulation, when the capsule shape is still a sphere and thus the stretching energy vanishes for the spherical shape. The spring constant k s is related to the 2D Young modulus via K 0 = 2k s / √ 3 [33]. The bending energy is included via a potential between normals on triangles sharing an edge where n α is the unit normal on triangle α, and α and β are triangles sharing an edge. The summation runs over all pairs of triangles. The bending rigidity κ is related to the coupling constantκ via κ = √ 3κ/2. This relation is obtained by calculating the bending energy of a regularly triangulated cylinder [33,37]. Networks with this type of stretching and bending potentials were studied in relation to the stability of membranes and icosahedral shells that contain defects [7,33] or defect scars [38]- [40] as well as for the deformation of icosahedral viruses [41]- [44] and the crumpling of elastic sheets [29]. Note that the spontaneous curvature is zero in this description such that the ground state is flat when only bending energy is considered.
In experimental systems the convenient method to pressurize the system is to apply an osmotic pressure. For instance, vesicles or capsules are prepared in some solution of macromolecules, after which they are put into a solution with a different osmolarity. Initially the pressure difference will be large, but by exchanging solvent (= volume) between bulk and the inside of the vesicle the osmotic pressure difference decreases. In this way neither the pressure nor the volume is constant in this system. Such a situation was considered by Linke et al, who studied the osmotically driven passage of a fluid vesicle through a narrow pore [45]. The volume after the pressure instability can be controlled by the amount of solute molecules inside the vesicle. For a very small amount of solute in the vesicle, no pressure will build up inside the vesicle, there will be an almost complete collapse and a constant-pressure ensemble is appropriate to describe the system. On the other hand, for a large amount of solute in the vesicle the concentration will rapidly increase as the volume becomes smaller and it is not difficult to adjust the volume to any desired value. We consider the two limiting cases. The first one is to add a pressure-volume term to the Hamiltonian so that the external pressure is the control parameter to drive the system through the collapse transition. The second one is to fix the volume of the vesicle or capsule and reduce it in a controlled sequence of steps, as is shown in figure 3. The volume change between subsequent steps (δV ) as well as the number of MC-steps per volume (N s ) determines the compression rate τ v = δV /N s V 0 . By the freedom of choosing δV and N s independently, there are still an infinite number of choices to establish a certain compression rate. The constant-volume constraint is implemented via a potential, where V is the volume of a specific configuration and V the desired volume. Note that for k v too large, the evolution of the system will be completely determined by this volume term. In the case that k v is too small, the volume term will not suffice to compress the system. Most simulations discussed here are done with shells of 5530 vertices and R 0 ≈ 20/r 0 , where r 0 is the average lattice constant in the triangulation. The temperature in the simulations is chosen to be so low that the effect of thermal fluctuations is negligible, i.e. the smallest bending rigidity employed is κ/k B T = 5800. For the volume constraint, k v r 6 0 /k B T = 100 is used.  We performed two types of simulations. In the first case, we prepared the initial state as a sphere with a single or with multiple indentations (two or six) of depth H symmetrically distributed over the sphere. After equilibration at fixed volume, that is when the total energy has reached a plateau value, we calculate the value of the total energy.
The natural way to link experimental systems to the simulations is via the Föppl-von Kármán number γ given in equation (1). Typical experimental values are R ≈ 3 µm and h = 10 nm such that γ ≈ 10 6 [12].

Data analysis
In the simulations, we have calculated the total elastic energy of the deformed sphere as a function of the volume. Moreover, we calculate the local mean curvature on each of the vertices using the algorithms detailed in [46,47]. From this we obtain the distribution function of the mean curvature.
Furthermore, we calculate the height-height correlation function [48] where h(r) = R(r i ) − R is the local deviation of the radius from the average radius R of the shell. We calculate h(r i ) and h(r j ) as well as the distance r i j between the points r i and r j along and projected on the average sphere as illustrated in figure 4. The data are collected in a histogram where the distance r is stored in discrete bins of width δ = 0.75r 0 and the value of h(r i )h(r j ) is averaged for each bin. The correlation function shows an oscillatory character where following [48] the first zero crossing λ characterizes the wavelength of the pattern. The roughness of the shell is expressed by the amplitude of the correlation function. pressures for which the shell has collapsed. The black line that separates the spherical and the collapsed state is to guide the eye. In all these simulations K 0 r 2 0 /κ = 10 and κ 0 /k B T = 5800 whereas 7.5 < R 0 /r 0 < 55.

Pressure-induced collapse
We first performed simulations in which a spherical shell is exposed to an external pressure by adding a pressure volume term to the Hamiltonian. Figure 5(A) shows the radius of the shell (after equilibration) as a function of the applied pressure. For small pressures the radius slowly decreases with increasing pressure. At the critical pressure, the shell collapses and the radius rapidly drops. The critical pressure, figure 5(B), nicely follows the expected 1/R 2 0 dependence predicted by equation (2). The prefactor of 3.55 that we find is slightly smaller than the theoretical value of 4. This might on the one hand be due to the subtleties involved in the conversion between the coupling constantsκ and κ (as pointed out above) and on the other hand be caused by soft spots, vertices that are slightly softer than average due to, for instance, a locally lower area density of vertices.

Regular deformations
We study the regular deformation of a spherical shell by making single or multiple indentations of the same depth. The problem of the elastic energy of a single indentation was already studied many years ago [32], and we will just give the result here. Consider a spherical shell with an indentation of depth H as indicated in figure 1. The volume change V is just twice the volume of the spherical cap,  (14).
After elastic relaxation, the rim of the inverted cap will be smoothly curved and most of the elastic energy is localized in this rim. Since the spontaneous curvature is zero, the inverted spherical cap of course has the same energy as the original cap. Note that for the case of an equilibrated shell, equation (9) might not be exactly applicable anymore since the indentation might lead to a global deformation of the shell. However, we assume that the deviations are small. Assuming that the rest of the sphere is not distorted, the elastic energy associated with an indentation of depth H is [32] where Q 0 is a numerical prefactor. In order to find the energy as a function of the relative volume change, we have to invert equation (9); this cannot be done exactly, but a good approximation can be obtained using a series expansion. For small indentations we find that where V = 2V c and V 0 = 4π R 3 0 /3 so that This is the energy of a single indentation. If we now distribute the volume that is taken out over N i smaller indentations, the total elastic energy of these N i caps becomes so that for small indentations the energy grows like N 1/4 i . For larger indentations, we employ a series expansion in where X = V /V 0 and the coefficients A 1 -A 6 are listed in table 1. Using the truncated expansion (equation (10)) for the energy, we can immediately write as the energy for N i indentations with total volume V . Here X = X/N i . Figure 6(A) shows H/2R 0 as a function of V /V 0 as calculated from equation (14). The agreement between equation (14) and the series expansion is very good, so that for our purposes the number of terms in the expansion is sufficient. Figure 6(B) shows the energy as a function of V for various numbers of indentations as calculated using equation (15). For all deformations the ground state is a single indentation. We note that this result is in contradiction with the results of Quilliet [49], who finds that for a large enough volume change (for V /V 0 0.1) multiple indentations become more favourable.
As an ultimate test, we calculate the energy of configurations with 1, 2 and 6 indentations for 0.1 < V /V 0 < 0.5 from simulations. In these simulations, the start configuration is created by pre-indenting the sphere with 1, 2 or 6 indentations such that the overall volume is the desired one. Subsequently, the configurations are allowed to relax in a Monte Carlo simulation keeping the volume fixed. In figure 7(A), the results are shown for a relatively small Föppl-von Kármán number γ = 1.6 × 10 3 . At first glance it is surprising that all data collapse onto a single curve. However, after a visual inspection of the configurations it turns out that after equilibration the systems with 2 and 6 indentations have in fact relaxed to configurations with a single indentation such that the data coincide. This type of relaxation of a state with six indentations towards a state with only one indentation is illustrated in Movie 1 for a shell with an even smaller Föppl-von Kármán number γ = 700. It is very clear that in the beginning one of the indentations starts to grow at the expense of another in an Ostwald ripening type of coarsening, and this process repeats itself until only one indentation is left.
In figure 7(B), the results for γ = 1.6 × 10 5 are shown. In this case, the data do not collapse and visual inspection of the configurations shows that for 2 and 6 indentations these are still present. Clearly, these two results demonstrate that for these reduced volumes the single indentation always has lower energy than multiple indentations. Moreover, the lines in figure 7(B) are obtained by using equation (15) with the same prefactor Q 0 = 24 for all three data sets. This prefactor is undetermined in the scaling argument and the reason that figures 7(B) and 8 do not have the same energy scale. When the data for a single indentation for the two values of γ are compared, we find that the curves are similar but not the same, indicating no large 'energetic' differences between an indentation for small and large γ . However, for small γ the scaled energy is systematically larger. We note that for small γ the shape of the shell is slightly distorted also in regions far away from the indentation so that the validity of equation (15) is questionable. At the same time, for large γ the indentation becomes polygonal, which may in turn lead to deviations from equation (15) as well. We will return to this point in the next section.

Shape of the indentation and ridge scaling
As observed in experiments [23] and simulations of adhesive capsules [50], indentations can become polygonal where the sides of the polygon resemble stretching ridges very similar to the ones that are found in crumpled elastic sheets. However, in theoretical calculations of the energy of an indentation, the possibility of the formation of such stretching ridges has not been taken into account. We present here simulation results for the energy of an indentation of a certain depth (or equivalently a shell with a certain volume) as a function of K 0 . The increase of K 0 leads to a proportional increase of γ , and for large enough γ the rim of the indentation is expected to change shape, i.e. the rim becomes polygonal. If ridge scaling is applicable, we expect the scaling of the energy to change from E ∼ κγ 1/4 for a smooth indentation (compare equation (10)) to E ∼ κγ 1/6 for ridge scaling [25,26,29]. In order to test this, we have calculated the total energy of a single indentation as a function of γ . In figure 8, we show the energy as a function of γ for four reduced volumes. In all cases we see, as expected, a gradual change of the scaling from E ∼ γ 1/4 for small γ to a smaller exponent; however the dependence is found to be much weaker than γ 1/6 . This is consistent with the energy scaling observed by  Lidmar et al for icosahedral shells [7]. In their calculations of the energy as a function of the Föppl-von Kármán number, very weak γ -dependence is observed for intermediate γ -values separating the low and high γ limiting scaling behaviour of the energy. Here, the system sizes are too small to verify the large-γ limit in detail. For all configurations we determined the structure of the rim (smooth or polygonal) and this information is included in figure 8 as well. A slightly more detailed analysis of our data shows the empirical scaling γ p 7000( V /V 0 ) −3/4 for the transition between a smooth and a polygonal rim. In figure 9, we show the scaling of the stretching energy and bending energy separately. In figure 9(A), we see that the stretching energy rapidly increases with increasing γ when no ridges are present; however, when ridges are formed, the stretching energy starts to decrease again. We note that for large γ the typical length scale d r 0 so that discretization errors are more likely to affect the properties of the shells. We therefore do not put too much value to the increase of the stretching energy for large γ as found in figure 9(A). The bending energy ( figure 9(B)) shows a monotonic increase with increasing γ .

Dependence of the shape on Föppl-von Kármán number and compression rate
After a volume change δV , when the pressure exceeds the critical pressure, a homogeneous shell becomes unstable to small wavelength indentations, that is indentations larger than the thickness of the shell h ∼ √ κ/K 0 . The undulation amplitude grows and the wavelength is expected to coarsen until the new volume is reached and the shape has relaxed to a local minimum of the elastic energy. However, a shell in a simulation or experiment is never completely homogeneous. Therefore, an indentation can already appear at soft spots slightly below the critical pressure.
The number of such nucleation sites will increase with faster compression rate. Furthermore, for slow compression rates, the system has time to relax into its local minimum, whereas for high compression rates, the shell will not reach the relaxed state before the next volume change δV occurs. This drives the system further out of equilibrium. This means that during the process of volume reduction there is a competition between relaxation, i.e. the number of indentations decreases, and a deepening of the existing indentations. We present simulation results for three different Föppl-von Kármán numbers γ = 8 × 10 3 , γ = 8 × 10 4 and γ = 8 × 10 5 , in which the volume is reduced in a stepwise fashion as discussed before. For each of these systems, we performed simulations at three compression rates (τ v = 1.5 × 10 −6 , 1.5 × 10 −7 and 1.5 × 10 −8 ) to a final volume of V /V 0 = 0.48. For the largest compression rate τ v = 1.5 × 10 −6 , we studied for each of the Föppl-von Kármán numbers how the energy and the number of indentations relax at fixed volume after reducing the volume to V /V 0 = 0.48. In figures 10, 11 and 12, snapshots are shown for the three γ values and the three compression rates.
In the top row of figure 10 the evolution is shown for the fastest compression rate at γ = 8 × 10 5 . Below V /V 0 = 0.93 the shell becomes unstable and the appearance of indentations is observed. On further decrease of the volume the shell crumples progressively. For the smaller compression rate τ v = 1.5 × 10 −7 (middle row), the typical length scale is larger, that is, there are fewer indentations and the number of indentations decreases steadily upon further compression, as can be seen from figure 13(A). Finally for the smallest compression rate τ v = 1.5 × 10 −8 (bottom row), several indentations nucleate rapidly but almost instantaneously coarsening takes place and only three indentations remain when the system has reached the final reduced volume V /V 0 = 0.48. In movies 2 and 3 (available from stacks.iop.org/NJP/13/045020/mmedia), the process of nucleation and coarsening is shown for γ = 8 × 10 5 and for systems with τ v = 1.5 × 10 −7 and τ v = 1.5 × 10 −8 . The movies show the elastic energy in a Mercator projection. In movie 2, which shows the evolution at τ v = 1.5 × 10 −7 , the initial instantaneous nucleation of very many small indentations is observed whereas in time (volume reduction) a fraction of these indentations disappear until approximately half of the indentations are left. In movie 3, which shows the evolution of the system at slowest compression τ v = 1.5 × 10 −8 , one can see the nucleation of nine indentations, of which three survive the process of coarsening. In both cases the mechanism for coarsening seems to be Ostwald ripening rather than the fact that indentations merge. A more detailed analysis of the number of indentations as a function of reduced volume and Föppl-von Kármán number will be given below, after the more general features of the snapshots are discussed. For γ = 8 × 10 4 the snapshots are shown in figure 11. The top row shows again the fastest compression rate. A comparison of these snapshots with the equivalent ones (those of the top row) in figure 10 shows that the smaller γ implies a considerable smoothing of the folds and a reduction of the number of indentations. In movie 4 we show the formation and relaxation of indentations for this compression rate. As before, with the smaller compression rate τ v = 1.5 × 10 −7 (middle The formation and relaxation of the systems compressed at τ v = 1.5 × 10 −7 and τ v = 1.5 × 10 −8 are also shown in movies 2 and 3, respectively (available from stacks.iop.org/NJP/13/045020/mmedia). row) the system again forms fewer folds and finally, for the system with smallest compression rate τ v = 1.5 × 10 −8 , two indentations are nucleated, of which one remains. For this single indentation the increasingly polygonal character of the rim with deepening indentation can be seen very clearly. Finally, the snapshots for γ = 8 × 10 3 are shown in figure 12. Following the trend, this system also forms structures that appear crumpled when compressed fast (top row) and more regular and smooth when compressed slower (middle and bottom rows). Compared with the systems with γ = 8 × 10 5 and γ = 8 × 10 4 , the folds are smoother. Interestingly, the simulation at smallest compression rate here does not produce a single but a double indentation.
In figure 13, the number of indentations is plotted as a function of reduced volume. Here we only give the data for the two slowest compression rates, since for the largest compression rate the fold patterns are too fuzzy so that the quality of the images in the Mercator projection does not allow for such an analysis.
For the larger compression rate τ v = 1.5 × 10 −7 shown in figure 13(A), the number of indentations is much larger (by a factor of 5) than for the smaller compression rate τ v = 1.5 × 10 −8 shown in figure 13(B). In all cases the absolute number of indentations decreases approximately by a factor of 2 over the whole range of volume change. Moreover, for the system in figure 13(A) with τ v = 1.5 × 10 −7 the number of indentations is much larger for the system with γ = 8 × 10 5 than for the systems with γ = 8 × 10 4 and γ = 8 × 10 3 , which behave very similarly.
In figure 14, the elastic energy, normalized by the bending rigidity, is shown as a function of reduced volume for three Föppl-von Kármán numbers and three compression rates. Also included are data points for a single indentation.
Not surprisingly, it is found that for all γ values the shell with a single indentation is the lowest in energy. Furthermore, when the system is further driven out of equilibrium by increasing the compression rate, the elastic energy increases as well. For γ = 8 × 10 4 and γ = 8 × 10 3 the compression curves overlap for small V /V 0 . For these small reduced volumes the shells are still spherical and have not buckled yet. The buckling point moves towards smaller reduced volume for smaller compression rates since a smaller compression rate allows a more effective stress condensation followed by buckling. After full compression to V /V 0 = 0.48, when the volume is kept fixed, the energy slowly relaxes as is shown in figure 15(A). The system with γ = 8 × 10 3 value relaxes faster than those with γ = 8 × 10 4 and γ = 8 × 10 5 . In fact, for the longest time scales achieved in the simulations, we find that the system with γ = 8 × 10 3  relaxes to a state with two indentations, γ = 8 × 10 4 to a state with three indentations and the system with γ = 8 × 10 5 is still completely crumpled. During relaxation, the systems coarsen, with small indentations disappearing at the expense of larger ones. However, in none of the systems studied here have we observed full relaxation towards the system with only one indentation. The relaxation of the number of indentations after the state V /V 0 = 0.48 is reached is shown in figure 15(B) for γ = 8 × 10 3 and γ = 8 × 10 4 . For these cases the relaxation of the number of indentations follows an effective power law N ind ∼ t −0.375 as a function of time, with an exponent that is independent of γ for these two γ values. Figure 16 shows snapshots during the same relaxation process as displayed in figure 15. The images of the system with γ = 8 × 10 5 confirm that this system is still very crumpled and far from the ground state, whereas systems with γ = 8 × 10 3 and γ = 8 × 10 4 have relaxed substantially.
Finally we calculated the height correlation function and the mean curvature for the three γ values. First we consider how the systems with various γ and compression rates develop during compression from the spherical to the compressed state with V /V 0 = 0.48, and subsequently we study the relaxation for fixed volume when the system has reached V /V 0 = 0.48. Figure 17(A) shows the height correlation function for γ = 8 × 10 4 and four reduced volumes for the fastest compression rate τ v = 1.5 × 10 −6 . It can be seen that upon reducing    the volume the amplitude of the correlation increases (the effective thickness of the crumpled shell increases) and the typical wavelength of the indentations increases. From these correlation functions, as well as the ones for γ = 8 × 10 3 and γ = 8 × 10 5 , we extracted the typical length λ associated with the first zero of G(r ). Figure 17(B) demonstrates that λ increases with decreasing volume, which implies that relaxation continues during volume reduction. For small volume reduction V /V 0 0.2, there is little difference in the length scales for different γ -values.
Relaxation proceeds when the simulations are continued with fixed volume V /V 0 = 0.48, as can be seen in figure 18; this relaxation process is very similar for γ = 8 × 10 3 and γ = 8 × 10 4 .
An example of the distribution function of the mean curvature, in this case for γ = 8 × 10 3 and the compression rate τ v = 1.5 × 10 −6 , is shown in figure 19. For configurations with indentations, these distribution functions are in general characterized by two peaks: one at negative curvatures indicating the part of the surface that is indented, and one at positive  curvatures indicating the rims of the indentations. In figure 19, we see that for V /V 0 = 0.96 the curvature distribution has one peak centred around the radius of the original sphere. For the larger compressions V /V 0 = 0.63 and 0.48, the shell is buckled and two peaks in the distribution are visible.
We extracted curvatures corresponding to the locations of the peaks for all three γ values and plotted them as a function of reduced volume (A) or MC-steps (B) in figure 20. For the compression simulations (figure 20(A)), the peak at positive curvature shifts towards larger values when the system is more strongly compressed, which points to the formation of increasingly sharper folds. In the subsequent relaxation simulation ( figure 20(B)), the curvature decreases again.
Indeed, equation (3), which predicts that the typical length scale (curvature in this case) is proportional to γ 1/4 , is only valid at equilibrium. The curvatures in the simulations of compression and of relaxation do not show a clear γ -dependence. In order to verify equation (3) we return to the result of figure 9, which suggests that the scaling with γ 1/4 depends on the indentation depth and is only visible for small indentations, and analyze the curvatures as  Figure 21 shows the scaling of the curvature of the rim as a function of γ for an indentation of V /V 0 = 0.90. Here we find excellent agreement with the expected scaling for γ -values up to γ = 1.6 × 10 4 . For larger γ , where the rim becomes polygonal, the peak indicating the rim curvature disappears from the distribution and transforms into a broad shoulder. Moreover, for these large γ -values noticeable change of the curvature distribution with increasing γ is found and thus the γ 1/4 scaling breaks down. When the energy compression curves for a shell with a single indentation are reanalyzed in more detail, we do not find a measurable change in the shape of the curve when the indentation transforms from smooth to polygonal at small indentation depths.

Nucleation of a single indentation
The compression rates employed in the previous section were not small enough to nucleate single indentations. In order to study the transition from the spherical state to a state with one configuration, simulations were performed at very small compression rates. When a sphere is slowly compressed by just reducing its radius, the stretching energy will increase quadratically, where N sp ≈ 3N is the number of springs. For small compressions, equation (16) can be approximated by The energy of a single small indentation (equation (12)) is Balancing the energies for the two limiting regimes yields a scaling V cr /V 0 ∼ γ −3/5 for the transition from a spherical to indented configuration, where V cr denotes the transition volume, a result that was also obtained in [24]. In figure 22 the energy is shown for the transition from a spherical to a single indentation configuration for three γ -values. All curves display two different regimes: for small compression (A) a quadratic dependence on the reduced volume and linear in k s consistent with equation (16), and at larger compression (B) a dependence following equation (18). We find an overshoot of the quadratic regime similar to a supercooling. The location of the 'jump', which is due to the energy barrier between the two states, is found at a larger volume than V cr . The location of this jump depends on the compression rate and shifts towards the true transition when the compression rate becomes smaller.

Summary
Our simulations have shown that the shape of a shell that is collapsed by an external pressure depends on the deformation rate as well as on the Föppl-von Kármán number. The ground state is a shell with a single indentation. The scaling of the energy E ∼ γ 1/4 as derived by Landau and Lifschitz [32] and used frequently only holds for small indentations, where the rim of the indentation is a smooth circle. For larger indentations the rim becomes polygonal. The critical volume for the transition from smooth to polygonal depends on the Föppl-von Kármán number and scales like γ p 7000( V /V 0 ) −3/4 . When a shell is deformed by applying an external pressure or by volume reduction, the deformation proceeds in parallel with relaxation of the elastic energy. The relaxation leads to a decrease in the number of indentations. The relaxation rate depends on the Föppl-von Kármán number as well as on the compression rate that is used to generate the initial state of relaxation.