Scale-invariance in miniature coarse-grained red blood cells by fluctuation analysis

Abstract

To accurately represent the morphological and elastic properties of a human red blood cell, Fu et al. [Fu et al., Lennard-Jones type pair-potential method for coarse-grained lipid bilayer membrane simulations in LAMMPS, 2017, 210, 193–203] recently developed a coarse-grained molecular dynamics model with particular detail in the membrane. However, such a model accrues an extremely high computational cost for whole-cell simulation when assuming an appropriate length scaling – that of the bilayer thickness. To date, the model has only simulated “miniature” cells in order to circumvent this, with the a priori assumption that these miniaturised cells correctly represent their full-sized counterparts. The present work assesses the validity of this approach, by testing the scale invariance of the model through simulating cells of various diameters; first qualitatively in their shape evolution, then quantitatively by measuring their bending rigidity through fluctuation analysis. Cells of diameter of at least 0.5 μm were able to form the characteristic biconcave shape of human red blood cells, though smaller cells instead equilibrated to bowl-shaped stomatocytes. Thermal fluctuation analysis showed the bending rigidity to be constant over all cell sizes tested, and consistent between measurements on the whole-cell and on a planar section of bilayer. This is as expected from the theory on both counts. Therefore, we confirm that the evaluated model is a good representation of a full-size RBC when the model diameter is ≥0.5 μm, in terms of the morphological and mechanical properties investigated.

---

1 Introduction

The red blood cell (RBC) is the simplest and most well researched blood-borne cell, making it an ideal candidate on which to develop techniques in whole-cell computational modelling.^1–3 It is a highly deformable, “rubbery” cell, able to recover its shape after squeezing through very narrow capillaries.⁴ The RBC is primarily comprised of a 2-component membrane, enclosing a cytoplasm fluid interior. The cell membrane is solely responsible for the elastic response of the cell due to the entirely viscous nature of the cytoplasm,⁵ with the cytoplasm instead responsible for the volume incompressibility of the cell.⁶ The membrane is composed of a lipid bilayer and distinct cytoskeleton network bound to its inner surface, connected by transmembrane proteins. As the membrane thickness is much lower than the diameter of the whole cell, it has three-dimensional structure describable by two-dimensional elastic parameters.⁷ The resistance to bending of the cell is then characterised by the bending rigidity B. The lipid bilayer is essentially a two-dimensional fluid-like structure embedded in three-dimensional space, resistant to bending, but unable to sustain in-plane shear stress due to its highly diffusive nature.⁸ Conversely, the cytoskeleton resists shear deformation but has negligible bending resistance. Therefore, the bending rigidity of the whole-cell is dominantly produced by the bilayer, and should be measurable solely from a plane of bilayer lipids.^8–10

Measuring the bending rigidities of mesoscopic vesicles such as the RBC has been of experimental interest for decades, prompting development of many competing techniques.^11–15 When observed in vitro, the RBC membrane is seen to flicker, now understood to be its experiencing stochastic thermal fluctuations.¹⁶ The first quantitative analysis of this phenomenon was performed by Brochard and Lennon in 1975.¹¹ Using phase-contrast microscopy, they extracted the bending rigidity from the power spectrum of fluctuations, pioneering the methodology of fluctuation analysis (spectroscopy). To improve upon the measurement accuracy, fluctuation analysis has been under continued development ever since, and successfully applied to a wide range of mesoscopic vesicles.^17–23 When performed experimentally, a camera captures consecutive snapshots of the cell membrane in the focal plane of the microscope (so-called “contours”), allowing the bending rigidity to be calculated from analysis of their thermal undulations over time. In 1995, Strey et al.²⁴ approximated B = 4 × 10⁻¹⁹ J from the excitations of just four orthogonal points on the RBC membrane. Faucon and Mitov^21,25 introduced consideration of the reduced membrane tension to account for conservation of the vesicle surface area and volume. Using their method, Rodríguez-García determined B = 2.7 ± 0.6 × 10⁻¹⁹ J from healthy RBCs, and Gracia et al.¹⁴B = 2.25 ± 0.8 × 10⁻¹⁹ J from lipids extracted from the RBC membrane.

Computational RBC models have been developed under many different formalisms. To date, the most popular have been those most computationally efficient, namely continuum methods such as the finite element method (FEM), and the aggressively coarse-grained (CG), dissipative-particle-dynamics (DPD).^2,3,26 However, with the ever increasing power of High Performance Computing (HPC) facilities, CG models of higher resolution and complexity have seen increasing popularity; notably the coarse-grained-molecular-dynamics (CGMD) RBC model developed successively by Drouffe et al.,²⁷ Yuan et al.²⁸ and Fu et al.²⁹ In this CGMD formalism, the molecules comprising the bilayer, cytoskeleton and transmembrane proteins are explicitly represented as CG particles that interact through Lennard-Jones-like potentials and Hookean bindings. To give meaningful separation between the two distinct membrane components, the length-scale of the system is defined by the thickness of the cell bilayer (5 nm). The elastic properties are dictated by the broader physics of the complete system and not explicitly specified in the interaction potentials. Conversely, comparative DPD and continuum RBC models require the a priori knowledge of elastic moduli, made explicit within their model functions.^30,31 Therefore, the CGMD model is unique in allowing evolution of these properties, crucially enabling the testing of how elastic parameters change under new stimuli.

However, CGMD cell models such as this are very computationally expensive, due to their high resolution requiring a very large number of particles to comprise a full-scale cell. Therefore, they have historically had simulations restricted to only small patches of membrane.^8,32 Only recently has a full-scale CGMD RBC been modelled in its entirety³¹. Despite being an implicit-fluid model with notable care taken to ensure computational efficiency, simulating the full-scale cell (consisting of 3.2 million particles) for 100 000 time-steps on 20 CPU cores still took on the order of a day. An alternative approach has been to simulate cells in “miniature”, as has been done to date with the whole-cell CGMD model of Fu et al.^29,33 To increase computational efficiency, cells are built from a far smaller number of particles. However, as the same high resolution length scaling is maintained, the physical cell size is significantly reduced. To our knowledge, such CGMD models of RBCs have only been validated qualitatively by studying their shape at rest and under flow conditions, with no quantitative verification of this “miniature cell” approach yet conducted.

1.1 Aims

The present work gives, to the best of our knowledge, the first quantitative assessment of the “miniature cell” approach to the CGMD RBC model of Fu et al.²⁹ We do this by measuring the bending rigidity B through fluctuation analysis, confirming its scale invariability in the whole-cell, and dominance by the lipid bilayer. Firstly, the simulation methodology is outlined, before attempting the evolution of initially spherical cells of various model diameters to a biconcave stable state. The evolved whole-cell systems then have many contours recorded over an order of ms, allowing determination of B from fluctuation analysis. These results are then assessed against mechanical-cell and RBC theory along with past in vitro studies, to discuss the invariance of aspects of the model to the physical cell diameter.

2 Materials and methods

2.1 RBC model

The model used in this work is that of Fu et al.,²⁹ who built upon the lipid membrane model of Yuan et al.²⁸ The model is that for a CGMD, 2-component RBC, with explicit representations of the lipid bilayer, cytoskeleton, and internal and external fluid particles (see Fig. 1A). The lipid bilayer is represented as a one-particle-thick monolayer of CG spherical particles, each representing a large number of constituent lipids. In determining a distance dependant function for lipid-lipid interactions in the bilayer, it is challenging to find a form that produces the correct diffusion of particles. It has been shown that the classical 12-6 Lennard-Jones (LJ) potential only produces two membrane phases, a solid at low temperatures and gas at high temperatures.³⁴ At small separations the inter-particle forces are too strong to permit particle diffusion, and at large separations too weak to keep particles bound together. To provide the intermediate fluid phase necessary to allow such behaviour, a two branch interparticle function can be adopted.²⁸

Fig. 1 (A) Pictorial representation of the RBC membrane components, showing the cytoskeleton network attached to the lipid bilayer. Actin junction complexes (actin protofilament and protein band-4.1) connect the spectrin tetramers. The cytoskeleton is tethered to the lipid bilayer via transmembrane proteins – immobile band-3 at the spectrin–ankyrin binding sites and glycophorin at the actin junctional complexes. (B) Graphic of our pre-evolved (spherical) RBC, with each CG particle type shown in colour: (red) lipids, (orange) trans-membrane proteins, (green) spectrin tetramers, (black) junction complexes, (yellow) ankyrin, (purple) internal fluid and (grey) external fluid. The graphic is split into two mirror halves, to make the distinct particle types clearer visually. The particles described in (B) are highlighted in (A) by opaque circles of corresponding colour.

This work adopts the lipid–lipid interaction potential of Yuan et al.,²⁸ hereafter referred to as the Yuan potential. The Yuan potential has been shown to represent well the mechanical properties of a RBC bilayer membrane, including a diffusive fluid phase, due to the separation of attractive and repulsive branches.²⁹ It has an orientational dependence which allows the complex lipid hydrophobicity to be represented, being essential for the self-assembly of the bilayer in an aqueous environment.^8,35 The membrane properties of spontaneous curvature c₀, bending rigidity B and diffusivity D are conveniently characterised by three Yuan model parameters, θ₀, μ_Y and ζ respectively. θ₀ signifies the most energetically favourable angular configuration between particles, with μ_Y weighting the energy penalty for deviation away from this. ζ controls the slope of the attractive branch of the potential. The potential also features the LJ-like parameters of length σ, energy well depth ε and cut-off radius r_c. See S1 of the ESI† for further detail on the formalism of the potential.

While bilayer–bilayer interactions are managed by the Yuan potential, all other particle–particle interactions operate through classic 12-6 LJ potential functions. Assuming a typical RBC curvature c₀ ∼ −0.5 m⁻¹, the membrane curvature is parameterised as sin θ₀ = −1.41 × 10^−3 36,37 (see S1 of the ESI† for detail on this calculation). A Hookean harmonic bond-potential is also used to couple the transmembrane proteins of the lipid bilayer to the cytoskeleton. The characteristic biconcave RBC shape can only be attained if the in-plane shear elastic energy is relaxed to zero, achieved in the biological cell by the constant structural remodelling of the spectrin network.⁸In silico this can be achieved by modifying the equilibrium bond length in the bond-potential so that the bond energy is initially zero, thus the cytoskeleton is initially stress free. The equilibrium bond length is made to be a variable corresponding to the initial bond lengths between each particle pair, rather than a constant as in the standard Hookean potential.²⁹ The cell is then able to relax to a biconcave shape. The model parameterisation used is summarised in Table 1, generally following Fu et al.²⁹ However, the particle masses set by Fu et al. fail to reproduce comparative elastic and diffusive properties in our simulations. Instead, better agreement is found with the mass of all membrane particles , and all other particles σ_m.

Table 1 Default pair potential parameters for each CG particle interaction type within the model, all given in LJ units. The mass of each CG particle type is also given. “Fluid” refers to both internal and external solvent. “Bilayer” refers to the lipids and transmembrane proteins. “Cytoskeleton” refers to spectrin, junction complex and ankyrin particles. In our implementation, some properties of the cell were found to be poorly represented by particular values from Fu et al.²⁹ Alternate values used here are those marked by a *

Particle interaction	Pair potential	σ	ε	r cut	ζ	μ Y	sin θ0
Fluid–fluid	LJ	2.7	0.2
Fluid–cytoskeleton	LJ	1.0	0.2
Cytoskeleton–cytoskeleton	LJ	1.0	0.2
Bilayer–bilayer	Yuan	1.0	1.0	2.6	4	3	−1.41 × 10^−3*
Bilayer–cytoskeleton	LJ	1.0	0.2
Bilayer–fluid	LJ	1.0	0.2

---

Particle type	m ^M
Bilayer, cytoskeleton	0.5*
fluid	1.0*

---

Simulations are run utilising the LAMMPS molecular dynamics coding package,³⁸ operated as a library within Python. LAMMPS handles the thermodynamic evolution of the system, while particular biophysical calculations are performed in the parent Python code. Nose–Hoover algorithms are used for thermostat and barrostot thermodynamics, as documented by Shinoda et al.³⁹ Initial particle configurations are input from an independent Python code which generates a 3D “supercell” volume containing the configuration of pseudo-particle types as required (see Fig. 1B), alongside particle classifications interpretable by LAMMPS. Simulations are performed in the system of non-dimensionalised LJ units. However, to compare results with past in vitro studies, quantities must then be converted to SI units. In the presented formalism, each variable i has an associated dimensional conversion parameter σ_i which relates the non-dimensional “model” LJ unit (denoted i^M) to “real” SI unit (denoted i^R). For example, the conversion of a length r from LJ units to meters is denoted r^R[m] = σ_r[m]r^M. All simulations were run on 1–4 nodes of the local supercomputer, with each node having two 14 core 2.4 GHz Intel E5-2680 v4 (Broadwell) CPUs, and 128 GB of RAM. See S2 and S3 of the ESI† for detail on the conversion of each physical unit and benchmark of the model respectively.

2.2 Planar fluctuation analysis

As RBC lipid membranes have a thickness much less than the diameter of the whole cell, they are classically describable as thermally fluctuating 2-dimensional sheets, with elastic energy dependant on their bending rigidity B and tension Σ.¹² Furthermore, the lipid bilayer dominates the bending rigidity of the whole cell, as the internal plasma is entirely viscous and the cytoskeleton has negligible rigidity compared to the bilayer. Therefore, the bending rigidity of the whole cell should be representable solely from a planar patch of lipid membrane.^7,8,32

The bending rigidity of a planar membrane patch can be calculated by planar fluctuation analysis, from its height fluctuation spectrum.¹² In the Monge representation, the height spectrum of the square planar membrane is given by h(r) ≡ h(x, y) = z_i − 〈h〉, being the offset in height from the average across the entire membrane 〈h〉.⁴⁰ The height of a cell z_i is defined simply by the mean height of all particles within its bounds. Performing a 2D fast Fourier transform (FFT) on h(r) then gives the fluctuation spectra

with power spectrum 〈|h_q|²〉 by the squared modulus of the resulting q-space height grid, averaged over many equilibrium configurations.⁴¹ In the continuum limit, this equates towhere q is the norm of the wavevector q = 2π(n_x, n_y)/L, and L is the side length of the patch. The bending rigidity of the bilayer can then be simply determined from a fit of the mean-squared-amplitudes against q on a tensionless membrane (Σ = 0).

2.3 Whole-cell fluctuation analysis

Once each cell has been equilibrated to its stable state, vesicle fluctuation analysis can be performed to determine the bending rigidity of the whole-cell. Simulations of the stable-state cells record the passive thermal fluctuations within series of contours: two-dimensional slices of membrane as seen from a top down view. As the biconcave regions were prompted to develop through the z axis, contours are taken in the (x, y) plane, thus having radial fluctuations measured in the plane of the disc. To remove stochastic rotation of the vesicle, the net angular momentum of membrane particles is set to zero at the start of each time-step. As fluctuation analysis is non-invasive, the orientation and centre-point of the cells then remain approximately fixed throughout the simulations. For each cell, contours are taken every 1000 steps (20 μs), for 1000 contours in total.

As we have direct access to the coordinate data of each CG particle constituting our simulated vesicle, a contour detection algorithm is not necessary. Instead, we determine a contour from the explicit position data of our CG lipids at a given time. Each lipid particle in the vesicle membrane at time t has Cartesian coordinates r(x, y, z, t). Converting to spherical coordinates in the equatorial cross-section (at θ = π/2) gives its radial position at angle ϕ = tan⁻¹(y/x) around the membrane disc, relative to the midpoint of the vesicle (x₀, y₀, z₀ = 0, t). The midpoint is simply defined as the centre of the maximum and minimum x and y coordinates of all lipids within the membrane, x₀ = (x_max − x_min)/2 and y₀ = (y_max − y_min)/2 respectively, approximately equivalent to the centre of mass.

The angular distribution of a contour 0 ≤ ϕ ≤ 2π is segmented into N_ϕ segments ϕ_i, each having radius ρ(ϕ_i, t) simply defined by the lipid within that segment at the largest radius. The ϕ-average radius of a particular contour is expressed as ρ(t) = 〈ρ(ϕ_i, t)〉. Fig. 2 shows an example polar plot of a contour. While our simulations are not subject to the resolution limit incurred through optical microscopy, a similar limit is instead imposed by our level of coarse-graining. As we maintain the same level of coarse-graining between all cells (defined by length unit σ_r = 5 nm) but simulate cells of different diameter, the larger cells comprise more particles and thus have a greater potential “resolution” in their contours. To maintain consistency between cells, we therefore make the number of radial arcs N_ϕ depend on D^M. As the contour circumference is πD^M and the lipid–lipid spacing is ∼σ_r, there will be approximately πD^M CG lipids cosmprising a contour. For the resolution to be sufficient, the number of contour segments must be a minimum of half the number of particles comprising that contour, meaning N_ϕ > πD^M/(2σ_r). Rounding the numbers slightly, we thus define N_ϕ = 3.6D^M.

Fig. 2 Example polar contour from a D^M = 100 cell, with radial distribution over N_ϕ = 3.6D^M = 360 angular segments. The inset plot shows the flattened distribution, with the ratio of the radius to the mean radius given against angle.

Once the simulations have completed and all the contours have been recorded, a separate Python code is used for their analysis. The method employed is that developed successively by Faucon et al.,²¹ Mitov et al.,²⁵ and Melerard et al.⁴² By decomposing the contour fluctuations into the spherical harmonics basis, their mean square amplitude can be derived to be

with reduced membrane tension .²¹ The mode n = 0 corresponds to variations in the mean vesicle radius, and n = 1 to variation in the centre of mass. Neither of these are relevant for the fluctuation analysis, so only modes n ≥ 2 are considered.

To establish a connection between the experimentally observable two-dimensional contour slices and three-dimensional model that leads to eqn (3), the angular auto-correlation function (ACF) is introduced. Use of the angular ACF acts to correct for this dimensional simplification, mitigating almost all associated error in the intermediate modes; the modes that are most useful to experimentalists.²¹ For a given contour, the angular ACF is given by⁴²

with angular segment about the contour γ, and reference radius R₀ = (3V₀/(4π))^1/3 taken as that defined by a perfect sphere of the same volume as the cell V₀. To a very good approximation, the cell volume and surface area remain constant.²¹ Expanding the ACF as a Legendre polynomial series, the coefficients B_n(t) can be related to 〈|U^m_n(t)|²〉 bywhere their mean is taken over many contour configurations. The constants B and can then be determined from a χ² fit of 〈B_n〉 against n. See S4 of the ESI† for a more thorough description of the analysis method.

3 Results

3.1 Cell equilibration

To test the miniature-cell methodology, whole-cells are generated at various model diameters D^M, with smaller cells comprising proportionally fewer CG particles. Before fluctuation analysis can be performed, the cells must first be evolved to their biconcave stable state. Each cell is initially generated as a spherical configuration of CG particles suspended within a fluid of water-like CG particles. The cell membrane is generated as two concentric spherical shells of bilayer and inner cytoskeleton, enclosing an internal CG fluid representing the cytoplasm. Largely following Fu et al.,²⁹ each particle type is then thermodynamically activated sequentially as follows:

• An isothermal–isobaric (NPT) ensemble is applied to the external water over 25 000 steps. Concurrently, the canonical ensemble (NVT) is applied to the internal fluid, only being under the influence of a thermostat.

• The spectrin, ankyrin and junctional complexes of the cytoskeleton are then equilibrated with the NPT ensemble over 25 000 steps.

• Finally, the lipids and transmembrane proteins in the bilayer are equilibrated using the NVT ensemble over 50 000 steps.

Each cell is equilibrated from its initially spherical state with time-step length τ = 0.02σ_t, pressure P^M = 0.05, and temperature ramping from T^M = 0.02–0.23.²⁹ Under the default parameterisation of Table 1, the shape evolution was found to introduce significant instabilities, and often result in the cells “popping”. To prevent this, some model parameters are changed specifically for these initial shape evolution simulations: the membrane particles have their mass increased to σ_m, and the energy well depth in the Yuan potential is increased to ε = 1.5, to temporarily strengthen the lipid bindings.

Membrane folding is then induced by reducing the initial number of internal fluid particles N_IF,0 to a final number N_IF (see Table 2). Small, equal and opposite forces are also applied to circular areas on each XY face of the cell to promote biconcave indents to manifest perpendicular to the XY-plane. Particles are deleted gradually to a final fraction of internal fluid particles n_IF = N_IF/N_IF,0, differing with cell size. The rate of compression has an effect on both the equilibrium shape and stability of the transition.^7,29 All cells have N_IF reduced at a constant rate of 3% every 5000 time-steps, determined to be most conducive to achieving a biconcave final state. The concave regions of the cells then develop gradually with the compression. Deviation from this rate results in alternative unwanted vesicle transitions such as to prolate, dumbbell rods, or inward or outward budding.⁷

Table 2 Equilibrated final states of each RBC cell size, with final ratio V/R³ achieved from chosen compression fraction n_IF. The number of bilayer to total particles is also given, as well as the number of steps run in the compression stage of the simulation

Final state
D ^M	50	75	100	125	150
V/R^3	2.02	1.90	2.08	2.03	2.00
n IF	0.68	0.49	0.42	0.37	0.32
Bilayer particles	8346	18 704	33 400	52 069	75 156
Total particles	34 944	114 592	269 296	521 437	896 703
Compression steps	55 000	85 000	95 000	105 000	115 000

---

Only cells with diameter D^M ≥ 100 (D^R ≥ 0.5 μm) were able to achieve a biconcave discocyte final state, with smaller cells instead relaxing to bowl-shaped stomatocytes (see Table 2). The degree of biconcavity in a healthy RBC can be characterised by the volume–radius ratio V/R³ = 1.57.³⁷ An optimal particle fraction n_IF for each cell is found by slowly deleting internal fluid particles until a target volume-radius ratio is reached. To achieve a consistent volume–radius ratio between cell sizes n_IF is found to be inversely proportional to D^M. Furthermore, below a critical ratio V/R³ ≲ 1.9 the internal fluid becomes unable to fill the region between the two enclosing membrane edges, closing the gap. This critical ratio does not appear to have a direct relationship with cell-size. However, the effect is more pronounced for the biconcave cells, being suppressed in the bowl-shaped D^M < 100 (D^R < 0.5 μm) cells. To maintain consistency between cells, n_IF is chosen such that cells are compressed to V/R³ = 2.0 ± 0.1, slightly higher than that of a physical RBC. However, cells of diameter D^M ≥ 100 still achieved morphologies closely resembling those of a healthy human RBC (see Table 2).

3.2 Planar fluctuation analysis

To measure the bending rigidity of the lipid bilayer specifically, planar fluctuation analysis is performed on a two-dimensional patch of CG lipids. A square patch of isolated lipid membrane of side lengths L^M = 120 is generated, suspended within a super-cell volume of dimensions {L_X, L_Y, L_Z}^M = {122.0, 122.0, 32.4}. Simulations of this system are run over 2 × 10⁶ steps of length τ = 0.001σ_t, equilibrated using a Berendsen barostat such as to achieve zero surface tension. To approximate the bending rigidity, the height fluctuation spectra of the patch is determined by eqn (1). On length-scales shorter than the thickness of the membrane, protrusion modes dominate the fluctuation spectra and eqn (2) breaks down.⁴⁰ To avoid this, the membrane is mapped to a discrete 2D grid of size N_q² such that the unit grid size l = L/N_q is greater than the thickness of the membrane: namely l = 2σ_r. Fig. 3 shows the resulting spectrum in the small wave-number region, where a fit of eqn (2) determines B_2D = 18.3k_BT, matching closely the result of B = 18k_BT by Fu et al.²⁹

Fig. 3 Plot showing the fluctuation spectrum of the planar lipid patch in the low q limit, with best fit giving the bending rigidity. Membrane configurations are taken every 1000 steps, and discarded if the surface tension is calculated to be greater than ±0.01σ_F/σ_r offset from zero.

A further use of this planar bilayer system is to determine the system time-scale, through measuring the mean-squared-displacement (MSD).⁸ Comparing the resolved diffusivity against a typical value for an RBC bilayer reveals the approximate simulation time scale of σ_t = 80 ns, matching closely with the previous value σ_t = 100 ns of Fu et al.²⁹ See S2.1 of the ESI† for detail on this calculation.

3.3 Whole-cell fluctuation analysis

To calculate the bending rigidity of the stable RBCs, the whole-cell fluctuation analysis is employed on the evolved cells. The coefficients 〈B_n(t)〉 are taken from the mean of 6 repeat simulations of each cell size, with standard error taken in that mean. In the bending dominated regime, the relaxation time of a mode n in a cell of radius R can be approximated aswith viscosity of the surrounding solvent η.^15,20 This shows the relaxation time of a mode to increase proportionally to the cube of the cell radius. Statistical significance in the lower modes thus becomes increasingly difficult to achieve at larger cell size. Therefore, we perform fits of eqn (5) across different modal regimes, with the minimum modes considered being higher with larger cell size (see Table 3). Fig. 4 shows the resulting spectra for the D^M = 50 (D^R = 0.25 μm) and 150 (0.75 μm) cells, and the bending rigidity as calculated by eqn (5). This clearly shows the stable modal region, gradually shifting to higher modes with increasing cell size. Modes n ≲ 10 have notably poorer statistics and are more dependant on the tension , causing an overestimation in B. There is also a notable cutoff at n∼ 30, beyond which B exponentially falls, indicating a resolution limit. The surface tension fit poorly, likely due to a combination of high error in the coefficients, and fits being performed at high modal regimes where the surface tension is insignificant. Therefore, we performed the fits fixing an assumed approximation .⁴³ To confirm that these modal ranges produced non-correlated coefficients 〈B_n〉, the time-lag ACF was generated over extended simulations of the D^M = 50 (D^R = 0.25 μm) and D^M = 150 (D^R = 0.75 μm) cells. Within these chosen spectral ranges, no notable correlation was seen within the first 1000 contours.

Table 3 Fitting results across all tested cell sizes within the given modal ranges

D ^M	D ^R (μm)	n min	n max	B (kBT)	Reduced-χ^2
50	0.25	4	24	17.0 ± 0.9	3.2
75	0.375	5	28	17.3 ± 0.8	8.6
100	0.5	5	29	18.9 ± 0.6	5.2
125	0.625	8	29	16.3 ± 0.7	1.6
150	0.75	10	29	18.0 ± 0.7	2.6
175	0.875	13	29	19.5 ± 1.0	1.8
200	1.0	14	29	19.6 ± 2.0	7.0

---

Fig. 4 Example results from the D^M = 50 (D^R = 0.25 μm) and D^M = 150 (D^R = 0.75 μm) cells, with (A) fitting over the spectrum of coefficients 〈B_n(t)〉, and (B) the bending rigidity B as calculated by eqn (5).

From Fig. 5, there does not appear to be any notable relationship between B and D within the bounds of error, indicating B to indeed be invariant to the physical cell diameter across the sizes tested. While there appears to be a slight increase at D^M > 150 (D^R > 0.75 μm), we do not believe this to be significant considering the much larger uncertainty in measurement of these largest cells. The mean rigidity across all cells is determined to be B_3D = 17.9 ± 0.4k_BT. However, there is more scatter about the mean than can be accounted for by the calculated error for multiple cell-sizes. This is most likely due to an under-representation of the error. For example, the complex evolved morphological differences between each cell-size may have introduced measurement error that has not been possible to account for. Given our energy unit σ_ε = 1.76 × 10⁻²⁰ J, our mean bending rigidity of B_3D = 17.9 ± 0.4k_BT = 3.15 × 10⁻¹⁹ J matches well with the literature on healthy human RBCs (2–9 × 10⁻¹⁹ J^14,16,44). Furthermore, the result from the analysis on the two-dimensional patch of bilayer (B_2D = 18.3k_BT) also falls within the margin of error of B_3D. Therefore, we have shown the model correctly reproduces the theory on both counts: (1) the bending rigidity is invariant to the cell diameter, and (2) the bending rigidity of the whole cell can be represented solely from a planar section of its bilayer.

Fig. 5 Bending rigidity B as calculated from the whole-cell fluctuation analysis at different cell diameter. The error is from the standard error in the mean from 6 repeats of each cell size. Blue arrows indicate the shape variation, with only D^M ≥ 100 (D^R ≥ 0.5 μm) cells being biconcave discocytes.

4 Discussion

The primary goal of this work has been to study the scale invariance of the employed CGMD RBC model, to help justify the use of miniaturised cells in future studies. As hoped, we have observed invariance to cell size in both the cell morphology and in the bending rigidity as calculated through fluctuation analysis. However, the CG nature of the model did introduce limitations on the morphology and contour resolution. Cells of different physical diameter were found to evolve to different total morphology (biconcave discocyte or stomatocyte), and cells of the same size would also develop varied non-axial deformities within these shapes. For example, cells would develop as non-axial discocytes where the cell thickness was non-uniform, or having concave regions of varying depth and shape. Only cells of D^M ≥ 100 evolved to the characteristic biconcave discocyte shape of a healthy RBC, with smaller cells relaxing to bowl-shaped stomatocytes (see Table 2). This variation in morphology is most likely due to the number of degrees of freedom available from the number of constituent CG particles. Smaller cells are comprised of fewer particles, thus will have fewer degrees of freedom for their configuration, making them less versatile to shape transitions. Therefore, we consider it reasonable to have found a lower D^M limit on the cells able to form the more complex biconcave shape, and a complex relationship between D^M and optimal internal fluid fraction n_IF.

Similarly, the CG nature of the model introduced a resolution limit on the fluctuation spectra. However, by applying the fluctuation analysis to contours of relative decomposition N_ϕ = 3.6D^M, a measurable range of modes was achieved for all cell sizes tested. Noise in the analysed spectra was notably higher for D^M > 150 (D^R > 0.75 μm), resulting in larger uncertainties in measurement of B. Considering the longer relaxation times in these larger cells this is not surprising, as they would require significantly longer simulation times to avoid producing poorer statistics. Computational limitations from poor model scaling (see S3 of the ESI†) made this unfeasible, thus also the reliable testing of cells of D^M > 200 (D^R > 1.0 μm). Considering this, the bending rigidity as measured by whole-cell fluctuation analysis is determined to be invariant to RBC scale based on the measured range of 0.25–1.0 μm diameter.

4.1 Experimental differences

When applying fluctuation analysis, there are some notable differences between the methodology when applied in silico to in vitro. Fundamentally, the bending rigidity of a vesicle is extracted by measuring the spatial undulations of its membrane over time. This is achieved experimentally through light microscopy, where a camera captures consecutive images of a vesicle over a chosen exposure time. As a microscope only allows visualisation in its focal plane, one is experimentally limited to capturing two-dimensional contour slices of the vesicle. The inherently three-dimensional analysis must then be simplified to the two-dimensional contours available, with the error this introduces being relative to the discrepancy of a vesicle from being spherical. As our simulations provide direct access to the fully three-dimensional position data of the constituent particles, this experimental restriction is not forced upon us. However, the purpose of this work is to verify that the utilised model correctly represents certain properties of the physical cell. By approaching this through comparison against past in vitro studies, it is best to emulate those experimental methods as closely as possible. Therefore, we still restrict ourselves to the same two-dimensional contour analysis, to best compare our results against such experimental studies.

There are many other benefits to not requiring a camera, notably in the potential resolution of a contour. As cameras have a finite resolution, there is then an upper limit to the fidelity at which a point on the membrane can be determined by optical microscopy. Experimentally, once a digital image of the cell has been captured, the radial profile of the membrane must be extracted by some form of “contour detection algorithm”. This is typically performed by classifying peaks in a grey-level intensity spectrum across the image. A fluctuation of wavelength shorter than the pixel size or resolution will then not be detectable. The early approaches of Faucon et al.²¹ and Duwe et al.⁴⁵ achieved accuracies of 100 nm and 250 nm respectively, of the order of the image pixel (being around 100 nm). By using more precisely the sigmoid-shape of the grey-level profile, Pecreaux et al.¹⁵ achieved an order of magnitude improvement, with radial measurements of the contour within 10 nm. Our having direct access to the particle positions bypasses these considerations. However, we have instead a comparable resolution-like limit from our level of coarse-graining. Each of our spherical CG lipids is σ_r = 5 nm in diameter, also with minimum inter-particle separation of σ_r. This then implies a comparable maximum resolution to that of Pecreaux et al.¹⁵ We observed a hard resolution limit at n = 30, lower than implied from our CG resolution. The wavelength of a mode is λ ∼ πD^M/n. The 30th mode thus has a wavelength λ ∼ 5σ_r for the D^M = 50 cell and λ ∼ 16σ_r for the D^M = 150 cell. If our resolution limit was σ_r = 5 nm as expected, we should theoretically be able to measure modes in their 100's. This discrepancy is most likely due to the discrete particle nature of the model.

Cameras introduce a further source of error in their aperture time (shutter speed) τ_SS: the length of time over which light is collected to take the image. The higher the mode of a fluctuation, the higher its frequency and the shorter its relaxation time. Therefore, τ_SS introduces a limit to the wave-numbers detectable, with only fluctuations having a longer lifetime being measurable (τ_n > τ_SS). From eqn (6), for a typical camera exposure time of τ_SS = 33 ms measuring fluctuations of an 8 μm diameter RBC of B = 10⁻¹⁹ J in a solvent of η = 10⁻³ kg m⁻¹ s⁻¹, this corresponds to a lower limit of n ∼ 4.¹⁵

While our simulations are not limited by a cameras shutter speed, they are notably more constrained by the total exposure time – the length of time over which contours are taken. The lower the mode, the longer its relaxation time thus the longer time it takes to gain sufficient statistics in its spectrum of configurations. While in vitro experiments are able to collect thousands of contours over minutes (2000 contours over 1 minute,¹⁵ 3–4000 contours over ∼2 minutes,²³ 15 000 contours over 10 minutes⁴²), we are heavily constrained by our model timescale σ_t. To maintain numerical stability in the simulations, the shortest time-step size possible was τ = 0.02σ_t = 1.6 ns. To simulate just 5 seconds of the D^M = 150 cell would then take up to 40 days to complete (see S3 of the ESI† for a benchmark). From eqn (6), it will take ∼1000 time-steps for the 5th mode in our D^M = 150 cell to complete a single relaxation. Therefore, it is not possible for us to gain comparable statistics in the lower modes.

Fundamentally, this combination of considerations of contour resolution, shutter speed and exposure time define the spectral-range of modes best considered in fluctuation analysis. Fluctuations of wavelength shorter than the optical resolution are not detectable, and those rapid relaxation times will be shorter than can be captured within a cameras shutter speed. In practice, this means that light microscopy is unable to reliably measure modes n ≳ 20.⁴⁶ Conversely, our comparatively high resolution and lack of a shutter speed allows reliable measurement at higher mode numbers (n ≲ 30). However, our very short exposure time heavily limits our ability to gain statistical significance in the lower modes, becoming even more pronounced with increasing cell size. Therefore, we analysed spectral ranges relative to the cell-diameter. Each cell had contours recorded every 1000 time-steps, for a total of 1000 contours. From eqn (6), this then equates to measuring at least 1500 full fluctuations in each assessed mode. However, even on a state of the art HPC, it was infeasible to simulate cells more than just an 1/8th the diameter of a full-scale RBC for a sufficient time to gain statistically meaningful results from this analysis. This highlights the computational challenges faced when attempting whole-cell CGMD simulation, and practical benefit of the miniature-cell approach.

5 Conclusion

The scale invariance of a CGMD model for a single RBC has been tested in silico, first qualitatively through the shape evolution, then quantitatively through fluctuation analysis. By evolving cells of various diameter from their initially spherical configurations, those of D^M ≥ 100 (D^R ≥ 0.5 μm) were found able to develop the characteristic biconcave discocyte shape of a healthy RBC. The passive thermal fluctuations of all evolved cells were then observed in their equatorial contours, and the bending rigidity calculated through fluctuation analysis using the angular ACF. The theory on RBC mechanics states that: (1) the bending rigidity of cells of micron-order diameter should be invariant to size; and (2) that the bending rigidity of the whole-cell should be equivalent to that of an isolated plane of its lipid bilayer. The model successfully represented the theory on both counts. No notable relationship was found between the rigidity and cell size over the range 50 ≤ D^M ≤ 200 (0.25 ≤ D^R ≤ 1.0 μm) tested, and the value B_2D = 18.3k_BT determined from the lipid plane fell within the margin of error of the mean from the whole-cell analysis B_3D = 17.9 ± 0.4k_BT. However, the very short exposure time computationally feasible at the larger cell sizes resulted in poor statistics for cells D^M > 150 (D^R > 0.75 μm), with significant increase in the measurement uncertainty. Therefore, we conclude the use of miniature CGMD cells of 100 ≤ D^M ≤ 150 best able to reproduce the morphological and elastic responses of an RBC in silico, in the properties tested. This finding supports the use of the miniature cell approach in further studies, with its considerable computational advantages opening up numerous possibilities in simulations of physically larger, more numerous and more complex CGMD cellular systems than have been performed to-date.

Author contributions

P. A. carried out all research, simulations, data analysis and writing of the article. S. H. and A. M. S. supervised the research and reviewed drafts of the manuscript, as the primary and secondary PhD supervisors.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank S.-P. Fu for their help in initial implementation of the model, and N. J. Brooks for helpful discussion about fluctuation analysis. This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol – http://www.bris.ac.uk/acrc/. P. A. was funded through an EPSRC studentship EP/R513179/1.

Notes and references

1. T. Ye, N. Phan-Thien and C. T. Lim, J. Biomech., 2016, 49, 2255–2266.
2. X. Li, P. M. Vlahovska and G. E. Karniadakis, Soft Matter, 2013, 9, 28–37.
3. X. Li, H. Lu and Z. Peng, Handbook of Materials Modeling: Applications: Current and Emerging Materials, 2020, pp. 2593–2609.
4. X. Li, Z. Peng, H. Lei, M. Dao and G. E. Karniadakis, Philos. Trans. R. Soc. A, 2014, 372, 20130389.
5. S. Henon, G. Lenormand, A. Richert and F. Gallet, Biophys. J., 1999, 76, 1145–1151.
6. C. T. Lim, M. Dao, S. Suresh, C. H. Sow and K. T. Chew, Acta Mater., 2004, 52, 1837–1845.
7. H. Yuan, C. Huang and S. Zhang, Soft Matter, 2010, 6, 4571–4579.
8. H. Li and G. Lykotrafitis, Biophys. J., 2012, 102, 75–84.
9. N. M. Geekiyanage, M. A. Balanant, E. Sauret, S. Saha, R. Flower, C. T. Lim and Y. Gu, PLoS One, 2019, 14, e0215447.
10. Z. Peng, X. Li, I. V. Pivkin, M. Dao, G. E. Karniadakis and S. Suresh, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 13356–13361.
11. F. Brochard and J. F. Lennon, J. Phys., 1975, 36, 1035–1047.
12. M. Mutz and W. Helfrich, J. Phys., 1990, 51, 991–1001.
13. E. A. Evans, Biophys. J., 1983, 43, 27–30.
14. R. Garcia, N. Bezlyepkina, R. Knorr, R. Lipowsky and R. Dimova, Soft Matter, 2010, 6, 1472–1482.
15. J. Pécréaux, H.-G. Döbereiner, J. Prost, J.-F. Joanny and P. Bassereau, Eur. Phys. J. E: Soft Matter Biol. Phys., 2004, 13, 277–290.
16. Y.-Z. Yoon, H. Hong, A. Brown, D. C. Kim, D. J. Kang, V. L. Lew and P. Cicuta, Biophys. J., 2009, 97, 1606–1615.
17. M. B. Schneider, J. T. Jenkins and W. W. Webb, J. Phys., 1984, 45, 1457–1472.
18. H. Engelhardt, H. P. Duwe and E. Sackmann, J. Phys. Lett., 1985, 46, 395–400.
19. I. Bivas, P. Hanusse, P. Bothorel, J. Lalanne and O. Aguerre-Chariol, J. Phys., 1987, 48, 855–867.
20. S. T. Milner and S. A. Safran, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 4371–4379.
21. J. Faucon, M. D. Mitov, P. Méléard, I. Bivas and P. Bothorel, J. Phys., 1989, 50, 2389–2414.
22. M. I. Angelova, S. Soléau, P. Méléard, F. Faucon and P. Bothorel, Trends in Colloid and Interface Science VI, Darmstadt, 1992, pp. 127–131.
23. J. Henriksen, A. C. Rowat and J. H. Ipsen, Eur. Biophys., 2004, 33, 732–741.
24. H. Strey, M. Peterson and E. Sackmann, Biophys. J., 1995, 69, 478–488.
25. M. D. Mitov, J. F. Faucon, P. Méléard, I. Bivas and P. Bothorel, in Advances in Supramolecular Chemistry Vol. 2, ed. G. W. Gokel, Jai Press, Greenwich, 1992, pp. 93–139.
26. V. Rajagopal, W. R. Holmes and P. V. S. Lee, Wiley Interdiscip. Rev.: Syst. Biol. Med., 2018, 10, e1407.
27. J. Drouffe, A. Maggs and S. Leibler, Science, 1991, 254, 1353–1356.
28. H. Yuan, C. Huang, J. Li, G. Lykotrafitis and S. Zhang, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 011905.
29. S.-P. Fu, Z. Peng, H. Yuan, R. Kfoury and Y.-N. Young, Comput. Phys. Commun., 2017, 210, 193–203.
30. D. A. Fedosov, B. Caswell and G. E. Karniadakis, Biophys. J., 2010, 98, 2215–2225.
31. Y.-H. Tang, L. Lu, H. Li, C. Evangelinos, L. Grinberg, V. Sachdeva and G. E. Karniadakis, Biophys. J., 2017, 112, 2030–2037.
32. H. Li and G. Lykotrafitis, Biophys. J., 2014, 107, 642–653.
33. M. Becton, R. D. Averett and X. Wang, Biomech. Model. Mechanobiol., 2019, 18, 425–433.
34. I. R. Cooke, K. Kremer and M. Deserno, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 011506.
35. H. Noguchi, J. Phys. Soc., 2009, 78, 041007.
36. C. Pozrikidis, Math. Med. Biol., 2005, 22, 34–52.
37. D. Hartmann, Biomech. Model. Mechanobiol., 2010, 9, 1–17.
38. A. P. Thompson, S. J. Plimpton and W. Mattson, J. Chem. Phys., 2009, 131, 154107.
39. W. Shinoda, M. Shiga and M. Mikami, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 134103.
40. Z.-J. Wang and D. Frenkel, J. Chem. Phys., 2005, 122, 234711.
41. M.-J. Huang, R. Kapral, A. S. Mikhailov and H.-Y. Chen, J. Chem. Phys., 2012, 137, 055101.
42. P. Méléard, T. Pott, H. Bouvrais and J. H. Ipsen, Eur. Phys. J. E: Soft Matter Biol. Phys., 2011, 34, 116.
43. F. Jähnig, Biophys. J., 1996, 71, 1348–1349.
44. J. Evans, W. Gratzer, N. Mohandas, K. Parker and J. Sleep, Biophys. J., 2008, 94, 4134–4144.
45. H. P. Duwe, J. Kaes and E. Sackmann, J. Phys., 1990, 51, 945–961.
46. Y. Z. Yoon, J. P. Hale, P. G. Petrov and P. Cicuta, J. Phys.: Condens. Matter, 2010, 22, 062101.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1sm01542g

---