Stochastic bond dynamics facilitates alignment of malaria parasite at erythrocyte membrane upon invasion
Abstract
Malaria parasites invade healthy red blood cells (RBCs) during the blood stage of the disease. Even though parasites initially adhere to RBCs with a random orientation, they need to align their apex toward the membrane in order to start the invasion process. Using hydrodynamic simulations of a RBC and parasite, where both interact through discrete stochastic bonds, we show that parasite alignment is governed by the combination of RBC membrane deformability and dynamics of adhesion bonds. The stochastic nature of bondbased interactions facilitates a diffusivelike reorientation of the parasite at the RBC membrane, while RBC deformation aids in the establishment of apexmembrane contact through partial parasite wrapping by the membrane. This bondbased model for parasite adhesion quantitatively captures alignment times measured experimentally and demonstrates that alignment times increase drastically with increasing rigidity of the RBC membrane. Our results suggest that the alignment process is mediated simply by passive parasite adhesion.
Introduction
Malaria is a dangerous mosquitoborne disease which kills nearly 0.5 million of people every year (World Health Organisation, 2018). It is caused by a protozoan parasite of the genus Plasmodium and proceeds in several stages (Miller et al., 2002; Cowman et al., 2012; White et al., 2014). After about 10 days from the initial infection through a mosquito bite, an infected liver releases a large number of merozoites, eggshaped parasites with a typical size of $12\mu m$ (Bannister et al., 1986; Dasgupta et al., 2014), into the blood stream. The blood stage of malaria infection is a clinically relevant stage, where merozoites invade healthy red blood cells (RBCs) and multiply inside by utilizing the RBC internal resources. This intraerythrocytic development is essential for merozoites to be hidden from the immune system and avoid clearance. After about 48 hours post RBC invasion, infected RBCs are ruptured and new merozoites are released into the blood stream to repeat this reproduction cycle. Thus, RBC invasion by merozoites is crucial not only for parasite survival, but also for further multiplication.
RBC invasion by merozoites is preceded by three key events: (i) initial attachment, (ii) reorientation or alignment of the parasite such that its apex is facing the RBC membrane, and (iii) formation of a tight junction (Koch and Baum, 2016). The apex contains all required machinery to invade RBCs after the tight junction is formed (Cowman and Crabb, 2006). At physiological hematocrit levels with a volume fraction of RBCs close to 40%, initial attachment of merozoites can be considered almost immediate after their egress from infected RBCs. However, the initial attachment has a random parasite orientation, which rarely provides direct alignment of the apex toward the membrane required to start the invasion. This implies that the parasite alignment is an extremely crucial step for successful invasion, which needs to be completed within a couple of minutes, as after this time period merozoites generally lose their ability to invade RBCs (Crick et al., 2014). To facilitate parasite alignment, merozoites contain a surface coat of proteins, mainly GPIanchored, which can bind to the RBC membrane (Bannister et al., 1986; Gilson et al., 2006; Beeson et al., 2016). However, one of the main difficulties in the investigation of RBCparasite interactions is that exact receptorligand bindings remain largely unknown. Electron microscopy images (Bannister et al., 1986) of merozoites adhered to a RBC suggest that along with short bonds of length $\simeq 20\phantom{\rule{thinmathspace}{0ex}}nm$, connecting the two cells, there exist much longer bonds of lengths up to $150\mathrm{nm}$, which may play an important role in early stages of merozoite adhesion to the RBC membrane. Furthermore, these long bonds have a much lower density than short bonds. Even though adhesion kinetics of such bonds remain unknown, recent optical tweezers experiments (Crick et al., 2014) indicate the adhesion force of spent merozoites to the RBC membrane to be within the range of 10 to 40pN.
Another important aspect during merozoite alignment is the deformation of the RBC membrane. Dynamic membrane deformations of various magnitudes are often observed (Dvorak et al., 1975; Gilson and Crabb, 2009; Glushakova et al., 2005; Crick et al., 2013) and are thought to aid in the alignment process (Weiss et al., 2015; Hillringhaus et al., 2019). Recent livecell imaging experiments show a positive correlation between RBC deformations and eventual merozoite alignment (Weiss et al., 2015). Most merozoites that successfully invade RBCs induce considerable membrane deformations, while the invasion success is much less frequent without preceding RBC deformations. Furthermore, these experiments lead to an estimate of an average alignment time of about $16\mathrm{s}$ (Weiss et al., 2015). A recent simulation study by Hillringhaus et al., 2019, with RBCparasite adhesion modeled by a homogeneous interaction potential, has confirmed the importance of membrane deformations, which facilitate parasite alignment through its partial wrapping by the membrane. However, this model shows static (not dynamic) membrane deformations and leads to average alignment times of less than $1\mathrm{s}$, indicating that an essential aspect of the alignment process has not been captured. Another speculation is that dynamic membrane deformations are induced actively by merozoites through changing locally the concentration of Ca+ ions (Lew and Tiffert, 2007; McCallumDeighton and Holder, 1992). This proposition has been confronted by recent experiments (Introini et al., 2018), which show that calcium release by parasite starts only at the invasion stage. Therefore, RBC membrane deformations are potentially induced by a passive mechanism, such as parasite adhesion.
In this paper, we focus on the passive compliance hypothesis (Introini et al., 2018) which assumes that RBC deformations and parasite alignment result from parasite adhesion interactions rather than from some active mechanism. Thus, our central question is whether parasite alignment can be explained purely by the passive compliance hypothesis. In contrast to the recent simulation study by Hillringhaus et al., 2019, where RBCparasite interactions are represented by a laterally smooth potential, the adhesion model presented here is based on discrete stochastic bonds between parasite and RBC membrane. This is a key step toward a realistic description of RBCmerozoite adhesion, since it eliminates the major shortcomings of the previous potentialbased model such as unrealistically fast alignment times and the absence of dynamic membrane deformations. Even though receptorligand interactions which determine parasite alignment are largely not known, our bondbased interaction model still incorporates a few experimental details such as the range of adhesion interactions and density of different agonists (Bannister et al., 1986). In particular, bonds of different lengths, that is long and short twostate bond interactions, are employed in the model. The bondbased parasite adhesion model generates an erratic motion of the parasite at the RBC membrane, visually similar to that observed experimentally (Weiss et al., 2015). Furthermore, it results in alignment times which agree quantitatively with those measured in experiments (Weiss et al., 2015; Yahata et al., 2012) and confirms the importance of membrane deformations for successful parasite alignment. The model is also used to investigate the effect of various adhesion parameters, such as bond extensional rigidities and kinetic rates, and ligand densities, on the parasite alignment process. Future investigations with this model can consider more realistic scenarios such as parasite adhesion and alignment under blood flow conditions.
The article is organized as follows. First, we introduce and calibrate our hydrodynamic model, where simulation parameters are tuned to quantitatively match several characteristics of the parasite motion at the RBC membrane from available experimental data by Weiss et al., 2015. Then, RBC membrane deformations and alignment times are investigated for this reference parameter set and several cases of altered bond kinetics and rigidities, and ligand densities. Finally, the effect of membrane stiffness on alignment times is studied.
Results
The RBC membrane is modeled as a network of ${N}_{\mathrm{rbc}}=3000$ vertices that are distributed uniformly on the membrane surface and connected by ${N}_{\mathrm{s}}$ springs (Gompper and Kroll, 2004; Fedosov et al., 2010a; Fedosov et al., 2010b; Fedosov et al., 2014). Our RBC membrane model incorporates elastic and bending resistance, and its biconcave shape is obtained by constraining the total surface area and enclosed volume of the membrane. Similar to the RBC, a parasite is modeled by ${N}_{\mathrm{para}}=1230$ vertices distributed homogeneously on its surface. The egglike shape of a merozoite (see Figure 1a) is approximated as (Dasgupta et al., 2014; Hillringhaus et al., 2019)
where $\mathrm{R}}_{\mathrm{a}}=1.5\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ and $\mathrm{R}}_{\mathrm{b}}=1.05\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ are diameters along the major and minor axes of the parasite, respectively. The parasite is much less deformable than the RBC, as no deformations of parasite body are visible in experiments (Weiss et al., 2015; Crick et al., 2014). Therefore, the merozoite is considered to be a rigid body, whose dynamics can be described by equations involving force and torque on the parasite’s center of mass and directional vector (Heard, 2006).
Both RBC and parasite are immersed in a fluid and the hydrodynamic interactions are modeled by the dissipative particle dynamics (DPD) method (Hoogerbrugge and Koelman, 1992; Español and Warren, 1995). The interaction of parasite and RBC membrane has two components. The first component corresponds to an excludedvolume repulsion to prevent an overlap between the two cells, which is modeled by the repulsive part of the LennardJones (LJ) potential with a minimum possible distance $\sigma =0.2\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. The distance $\sigma $ can be considered as an effective membrane thickness of a surface constructed from overlapping spheres with a diameter $\sigma $. Generally, $\sigma $ depends on the resolution length of both the RBC membrane and parasite (about $0.2\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ in our models) and is chosen large enough to guarantee no artificial membrane intersection or overlap between the cells. The effect of the precise value of $\sigma $ on simulation results is expected to be small and will be discussed later. The second interaction component represents adhesion which is modeled by discrete dynamic bonds between RBC and parasite vertices. Each parasite vertex represents one of the two different types of ligands: (i) long ligands with an effective binding range $\ell}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}=100\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m$ and (ii) short ligands with an effective binding range $\ell}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}}=20\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m$. Both ligand types are distributed randomly at the parasite surface with fixed ligand densities ${\rho}_{\mathrm{long}}$ and ${\rho}_{\mathrm{short}}$, such that their sum ${\rho}_{\mathrm{long}}+{\rho}_{\mathrm{short}}$ is equal to the parasite vertex density ${\rho}_{\mathrm{para}}$. Receptors for ligand binding are modeled by RBC vertices, each of which can bind only a single ligand, irrespective of its type. Due to the effective membrane thickness characterized by $\sigma $, long and short bonds can be formed by bound long and short ligands if the distance between RBC and parasite vertices is smaller than ${\mathrm{\ell}}_{0}+{\mathrm{\ell}}_{\mathrm{eff}}^{\mathrm{long}}$ and ${\mathrm{\ell}}_{0}+{\mathrm{\ell}}_{\mathrm{eff}}^{\mathrm{short}}$, respectively, where ${\ell}_{0}={2}^{1/6}\sigma$ is the equilibrium spring length that corresponds to the cutoff of repulsive interactions. Note that existing bonds are allowed to stretch beyond their effective binding ranges, see section ‘Methods and models’ for more details.
To relate simulation units to physical units, a basic length scale is defined as the effective RBC diameter ${D}_{0}=\sqrt{{A}_{0}/\pi}$ (${A}_{0}$ is the membrane area), an energy scale as ${k}_{\mathrm{B}}T$, and a time scale as RBC membrane relaxation time $\tau =\eta {D}_{0}^{3}/\kappa $, where $\eta $ is the fluid viscosity and $\kappa $ is the bending rigidity of the membrane. All simulation parameters in model and physical units are given in Tables 1 and 2. Average properties of a healthy RBC correspond to $\mathrm{D}}_{0}\simeq 6.5\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ with $\mathrm{A}}_{0}=133.5\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}{}^{2$ (Evans and Skalak, 1980) and $\tau \approx 0.92\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}$ for $\kappa =3\times {10}^{19}\phantom{\rule{thinmathspace}{0ex}}\mathrm{J}$ (Evans, 1983; Fedosov et al., 2010a) and $\eta =1\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{P}\mathrm{a}\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}$.
To better understand the effect of various adhesion properties on parasite alignment, several parameters such as bond formation and rupture rates, bond rigidity, and ligand densities are varied. For each fixed parameter set, a number of simulations are performed and the results are combined and/or averaged, which is necessary due to the stochastic nature of bondbased interaction as well as thermal fluctuation effects within the fluid. Note that each simulation is performed for a different random choice of parasite vertices which represent long and short ligands, while their densities remain fixed, see section ‘Methods and models’.
Calibration of RBCparasite interactions
A parasite adhered to the RBC membrane exhibits visually an irregular diffusivelike motion observed experimentally (Weiss et al., 2015), which is controlled by the ligand densities ${\rho}_{\mathrm{long}}$ and ${\rho}_{\mathrm{short}}$, bond rigidities ${\lambda}_{\mathrm{long}}$ and ${\lambda}_{\mathrm{short}}$, and the bond formation (${k}_{\mathrm{on}}^{\mathrm{long}}$, ${k}_{\mathrm{on}}^{\mathrm{short}}$) and rupture (${k}_{\mathrm{off}}$) rates that are currently not known. Nevertheless, available experiments (Bannister et al., 1986) suggest that the number of short bonds in RBCmerozoite interaction is lager than the number of long bonds, which is reflected in the ligand densities ${\rho}_{\mathrm{long}}$ and ${\rho}_{\mathrm{short}}$ assumed for our parasite model (see Table 2). To calibrate RBCparasite interactions, parasite dynamics at the RBC membrane (see Video 1) is quantified by its fixedtime displacement, which is measured by tracking the distance $\mathrm{\Delta}d$ traveled by the parasite at fixed intervals of time $\mathrm{\Delta}t$, see Figure 2a. Particle tracking is employed to measure $\mathrm{\Delta}d$ from available experiments (Weiss et al., 2015), where $\mathrm{\Delta}t$ is selected to be $1\mathrm{s}$, which is the time resolution of the experimental videos. Only time ranges, within which parasites remain visible and the RBC is not moving much, are included in the analysis.
Figure 2b compares experimental and simulated characteristics of fixedtime displacements for the interaction parameters given in Table 2. This set of parameters (further referred to as reference case) is obtained by varying ${\rho}_{\mathrm{long}}$, ${\rho}_{\mathrm{short}}$, ${\lambda}_{\mathrm{long}}$, ${\lambda}_{\mathrm{short}}$, ${k}_{\mathrm{on}}^{\mathrm{short}}$, ${k}_{\mathrm{on}}^{\mathrm{long}}$, and ${k}_{\mathrm{off}}$ until a good agreement between experimental and simulated parasite displacements is reached. However, the effective binding ranges of long and short ligands remain fixed at $\ell}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}=100\phantom{\rule{thinmathspace}{0ex}}\text{nm$ and $\ell}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}}=20\phantom{\rule{thinmathspace}{0ex}}\text{nm$ in this calibration procedure. The variance of experimental displacements in Figure 2b is larger than that in simulations due to a limited sample size of experimental data (20 samples). Note that this set of parameters is likely not unique, and other combinations of the parameters, which result in statistically similar parasitedisplacement characteristics, can probably be found.
To further characterize the parasite motion on the RBC membrane, the meansquared displacement (MSD) of the parasite’s center of mass is computed in simulations and shown in Figure 2c. At long enough times $t\gtrsim 3$ s, the parasite exhibits diffusivelike motion, indicated by a linear increase of the MSD curve with time. For shorter timescales, the MSD of parasite motion shows a transient anomalous subdiffusion, which may occur, for instance, in the case of sticky particle dynamics with alterations between sticking (i.e., stopping its motion for some time) and diffusing states (Saxton, 2007; Höfling and Franosch, 2013). The transient sticky dynamics is an appropriate description for an adhered parasite, where sticking periods correspond to time intervals within which no bonds are formed or ruptured. The diffusivelike dynamics is governed by the number of bonds ${n}_{b}$ and their on and offrates, as an adhered particle becomes slower and eventually gets arrested when ${n}_{b}$ is increased and the rates are decreased (Jana and Mognetti, 2019).
Parasite alignment
Recent experiments suggest that a successful RBC invasion strongly correlates not only with the distance between parasite apex and RBC membrane, but also with a perpendicular alignment of the merozoite toward the cell membrane (Koch and Baum, 2016). Furthermore, the junctional (invasion initiating) interaction range ${r}_{\mathrm{junc}}$ of the parasite’s apex is known to be around $10\mathrm{nm}$ (Bannister et al., 1986). Based on these observations, we define two quantities, (i) the apex distance ${d}_{\mathrm{apex}}$ from the RBC membrane, and (ii) the alignment angle $\theta $ that characterizes parasite orientation, both sketched in Figure 3a. Here, ${d}_{\mathrm{apex}}$ is defined as the distance between the parasite apex and the nearest membrane vertex,
the alignment angle $\theta $ as the angle between the parasite’s directional vector $\mathbf{\mathbf{n}}$ and the normal ${\mathbf{\mathbf{n}}}^{\mathrm{face}}$ of a triangular face whose center is closest to the apex,
Figure 3b,c shows distributions of apex distance ${d}_{\mathrm{apex}}$ and alignment angle $\theta $ for the calibrated RBCparasite interactions. Both characteristics are represented by distributions as the merozoite is very dynamic at the membrane surface. Minimum values of ${d}_{\mathrm{apex}}$ in Figure 3b correspond to the parasite’s apex being very close to the membrane (i.e., ${d}_{\mathrm{apex}}\approx \sigma $), whereas maximum values generally represent a configuration where the parasite is adhered sideways to the RBC. Furthermore, low values of $\theta $ in Figure 3c characterize the sideways adhesion orientation, while large values of $\theta $ represent a good alignment configuration. Note that an ideal merozoite alignment would be achieved if ${d}_{\mathrm{apex}}$ is less than $\sigma +{r}_{\mathrm{junc}}$ ($\mathrm{r}}_{\mathrm{j}\mathrm{u}\mathrm{n}\mathrm{c}}=10\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m$) and the alignment angle is $\theta \approx \pi $. Due to a discrete representation of the membrane, perfect alignment is unlikely, which requires to slightly relax these conditions. Therefore, we define a successful parasite alignment by the criteria
The choice of $0.8\pi $ in Equation 4 is also partially driven by the RBC discretization length of about $0.2\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. Half circumference of the parasite corresponds to $\pi {R}_{a}/2=2.36\phantom{\rule{thinmathspace}{0ex}}\mu \text{m}$, which is about twelve RBC discretization lengths. This means that our resolution in determining angle $\theta $ is close to $0.1\pi $, so that the window of $0.2\pi $ in the alignment criteria is large enough to avoid strong discretization effects.
In experiments, merozoite alignment times are measured as time intervals between initial parasite adhesion and the beginning of invasion (Weiss et al., 2015). Similarly, alignment time in simulations is calculated as the time required for the parasite to meet the alignment criteria in Equation 4 starting from an initial adhesion contact (i.e., formation of a few bonds). Figure 4b presents a distribution of alignment times from 86 statistically independent DPD simulations for the reference RBCparasite interactions in Table 2. The alignment times range between $1\mathrm{s}$ and $26\mathrm{s}$ with an average value of $9.53\mathrm{s}$. For comparison, the average alignment time was reported to be $16\mathrm{s}$ by Weiss et al., 2015, and the range of alignment times between $7\mathrm{s}$ and $44\mathrm{s}$ was found by Yahata et al., 2012, which agree reasonably well with our model predictions. Differences in alignment times between simulations and experiments are possibly due to a limited experimental statistics (e.g. only 10 samples in the study by Yahata et al., 2012) and/or selected model parameters, as the distribution of alignment times in our model can be altered by changing RBCparasite interactions. Therefore, further experiments and possible model improvements are needed to clarify the source of existing differences.
Note that the sample size (about 100) in simulations is limited by the computational cost. A single simulation, corresponding to a total physical time of about $26\mathrm{s}$, requires approximately 168 core hours on the supercomputer JURECA (Jülich Supercomputing Centre, 2018) at Forschungszentrum Jülich. Therefore, a direct bruteforce approach for the investigation of the effect of various parameters on the parasite alignment time is not feasible. To overcome this problem, MonteCarlo (MC) sampling (see section ‘Methods and models’ for details), which is based on a twodimensional probability map of parasite alignment characteristics (${d}_{\mathrm{apex}}$, $\theta $) illustrated in Figure 4a, is employed to determine the differences in alignment times for various parameter sets. Such a probability map is computed from several direct DPD simulations of RBCparasite adhesive interactions. Then, the MC procedure is used to model stochastic jumps between neighboring alignment states (${d}_{\mathrm{apex}}^{i}$, ${\theta}^{j}$) within the probability map, starting from a randomly selected initial state and continuing until the alignment criteria in Equation 4 are met, and the number of MC steps represents the alignment time. Distribution of alignment times ${t}_{\mathrm{n}}$ from the MC sampling is shown in Figure 4c for the reference parameter set. Clearly, the distributions obtained by direct (Figure 4b) and MC (Figure 4c) simulations are very similar, verifying the reliability of the MC approach. Note that alignment times ${t}_{\mathrm{n}}$ from MC sampling are measured in terms of MC steps, since MC simulations do not have an intrinsic timescale. The average alignment time for the reference parameter set is denoted as $\u27e8{t}_{n,\mathrm{ref}}\u27e9$ and assumed to be equivalent to $9.53\mathrm{s}$, the average alignment time from direct DPD simulations of RBCparasite adhesion. This implies that 10^{4} MC steps correspond to about $15\mathrm{s}$.
Membrane deformation and parasite dynamics
A recent simulation study by Hillringhaus et al., 2019 with a laterally homogeneous adhesion potential has demonstrated that the deformation of RBC membrane is crucial for a successful parasite alignment. Further, we show that ligand density, bond rigidity and kinetics not only control the parasite motion at the membrane surface, but also directly affect membrane deformation. To quantify the strength of membrane deformations, a change in total energy between the deformed state and the equilibrium state of the RBC membrane is computed as (Hillringhaus et al., 2019)
Figure 5 shows temporal changes in deformation energy, number of bonds, head distance, and alignment angle for the reference case. Two major contributions to the deformation energy (i.e. elastic stretching $\mathrm{\Delta}{E}_{\mathrm{sp}}$ and bending $\mathrm{\Delta}{E}_{\mathrm{bend}}$ energies) indicate that membrane deformation is very dynamic and has a strong variability in its intensity. This is due to the dynamic formation and dissociation of long and short bonds between the merozoite and RBC membrane.
An interesting observation is that the head distance and alignment angle in Figure 5 fluctuate around some average values, indicating that the parasite has a preferred orientation, which is consistent with a peak in the probability map in Figure 4a. To assess whether the most likely values of ${d}_{\mathrm{apex}}$ and $\theta $ are mainly determined by the egglike parasite shape, or also depend on the mechanical properties of the membrane, ${d}_{\mathrm{apex}}$ and $\theta $ distributions in Figure 3 for a deformable RBC are compared with those for the parasite adhered to a rigidified membrane (see section ‘Effect of RBC rigidity’) in Figure 3—figure supplement 1. Clearly, in the case of a rigid membrane, the preferred ${d}_{\mathrm{apex}}$ and $\theta $ values are determined by the egglike parasite shape, corresponding to a configuration with maximum adhesion area. In comparison to the deformable membrane (Figure 3), the peak in ${d}_{\mathrm{apex}}$ for the rigid RBC (see Figure 3—figure supplement 1) is shifted further away from zero. This indicates that the degree of wrapping has a significant effect on the preferred values of ${d}_{\mathrm{apex}}$ and $\theta $. Therefore, in addition to the egglike parasite shape, RBC membrane properties, such as bending rigidity, shear elasticity, and local curvature, affect the most probable values of ${d}_{\mathrm{apex}}$ and $\theta $. Furthermore, the fluctuations of ${d}_{\mathrm{apex}}$ and $\theta $ from their average values in Figure 5 represent parasite motion toward its apex or bottom due to stochastic bond dynamics. Thus, the parasite dynamics at the membrane can be described as a superposition of the rolling motion around its directional vector with a preferred orientation and intermediate fluctuations of parasite orientation toward its apex or the bottom. The rotational motion around the directional vector is preferred because it is not associated with a significant energy cost, while fluctuations in the orientation toward the merozoite’s apex or bottom have an energy penalty.
A further noteworthy result from simulations is that a successful alignment occurs more frequently in the concave areas of RBC dimples than at the convex rim of the membrane. This is due to the fact that the cell dimples have a favorable local curvature or a lower energy penalty for membrane wrapping (AgudoCanalejo and Lipowsky, 2015; Yu et al., 2018), which leads to a stronger parasite wrapping by the membrane, and thus a larger probability for successful alignment. Figure 5—figure supplement 1 shows that the merozoite forms more bonds in the dimples than at the RBC rim, confirming the positiondependent differences in membrane wrapping. Furthermore, our simulations show that merozoites move frequently into the dimple areas, starting from the initial rim contact, and remain there for the majority of simulation time. This behavior is again due to a more energetically favorable adhesion position within RBC dimples in comparison to the RBC rim. Energetically favorable parasite wrapping within the RBC dimples might be also advantageous for the subsequent entry into the cell.
The dynamic adhesive behavior of the parasite in the current stochastic bondbased model is in striking contrast to the previous adhesion model (Hillringhaus et al., 2019) based on a homogeneous interaction potential between the two cells, where no dynamic deformations were observed. A qualitative correspondence between these two models can be understood by considering a ratio ${k}_{\mathrm{on}}/{k}_{\mathrm{off}}=\mathrm{exp}(\mathrm{\Delta}{U}_{b}/{k}_{\mathrm{B}}T)$, where $\mathrm{\Delta}{U}_{b}$ is the binding energy of a single bond (Bell, 1978; Schwarz and Safran, 2013). Thus, the ratio ${k}_{\mathrm{on}}/{k}_{\mathrm{off}}$ directly controls the average number of bonds $\u27e8{n}_{b}\u27e9$ and the strength of adhesion (see section ‘Effect of bond properties on parasite alignment’), which are correlated with RBC deformation energy $\mathrm{\Delta}{E}_{\mathrm{rbc}}$. Similarly, in the parasite adhesion model with a homogeneous interaction potential (Hillringhaus et al., 2019), the strength of adhesion potential controls membrane deformations. Even though average membrane deformations can be compared for these two models, the stochastic bondbased adhesion model results in a very different diffusivelike dynamics of the parasite, which is governed by ${n}_{b}$ and the offrate ${k}_{\mathrm{off}}$ (Jana and Mognetti, 2019). A significant increase of ${n}_{b}$ and/or a decrease of ${k}_{\mathrm{off}}$ would lead eventually to parasite arrest (see section ‘Effect of bond properties on parasite alignment’), which can be compared well with the model based on a homogeneous interaction potential (Hillringhaus et al., 2019).
There exist three different timescales which might be relevant for the parasite alignment: (i) bond lifetime ${\tau}_{b}\simeq 1/{k}_{\mathrm{off}}$, (ii) membrane deformation time on the scale of parasite size ${\tau}_{p}\simeq \eta {R}_{a}^{3}/\kappa $, and (iii) rotational diffusion time of the parasite ${\tau}_{r}\simeq 8\pi \eta {R}_{a}^{3}/{k}_{\mathrm{B}}T$. These characteristic times are $\tau}_{b}\approx 0.013\phantom{\rule{thinmathspace}{0ex}}\text{s$, $\tau}_{p}\approx 0.011\phantom{\rule{thinmathspace}{0ex}}\text{s$, and $\tau}_{r}\approx 20\phantom{\rule{thinmathspace}{0ex}}\text{s$ computed from the model parameters given in Tables 1 and 2. There is a clear separation of timescales between ${\tau}_{r}$ and both ${\tau}_{b}$ and ${\tau}_{p}$, indicating that the rotational diffusion of the parasite is too slow to have a significant effect on merozoite alignment. Furthermore, ${\tau}_{b}$ and ${\tau}_{p}$ are comparable in magnitude, suggesting that both bond dynamics and membrane deformations are important for the alignment process. It is also interesting to note that the ratio ${\tau}_{p}/{\tau}_{r}={k}_{\mathrm{B}}T/(8\pi \kappa )\approx 6\times {10}^{4}$ depends only on the bending rigidity $\kappa $. This means that membrane deformation will always represent a dominating timescale over the rotational diffusion of the parasite, independently of the parasite size and the viscosity of suspending medium.
After the detailed analysis of parasite alignment, let us consider a possible influence of the effective membrane thickness, characterized by $\sigma $, on merozoite alignment. Figure 5—figure supplement 2 presents various alignment characteristics for $\sigma =0.15\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ and $\sigma =0.3\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ in comparison with the original choice of $\sigma =0.2\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. The simulation results indicate that the $\sigma $ value may affect the number of bonds between the RBC and parasite, and thus the degree of membrane wrapping. This result is not entirely surprising, as $\sigma $ also affects the binding range defined as $2}^{1/6}\sigma +{\ell}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}$ and $2}^{1/6}\sigma +{\ell}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}$ for long and short ligands, respectively. However, differences in alignment results are rather small for $\sigma =0.15\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ and $\sigma =0.2\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$, indicating that the choice for small enough $\sigma $ we made is appropriate. The case with $\sigma =0.3\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ exhibits a larger number of bonds and stronger membrane deformations than for $\sigma =0.2\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. Finally, note that fixedtime displacement characteristics of the parasite in Figure 5—figure supplement 2 remain nearly unaffected by the $\sigma $ value, because dynamical properties of the merozoite are mainly determined by the bond offrate, see the next section.
Effect of bond properties on parasite alignment
To better understand the dependence of merozoite alignment on bond kinetics, the offrate ${k}_{\mathrm{off}}$ is varied for two ratios ${k}_{\mathrm{on}}^{\mathrm{short}}/{k}_{\mathrm{on}}^{\mathrm{long}}$ of short and long bond onrates. Figure 6 presents the parasite’s fixedtime displacement, deformation energy, and average alignment times as a function of ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}$. A lower ratio of ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}$ (i.e. a lower ${k}_{\mathrm{off}}$) leads to stronger adhesion and thereby stronger membrane deformations (see Figure 6b and Video 2), consistently with the number of bonds shown in Figure 6—figure supplement 1. For small ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}$ values, membrane deformation energies can reach up to $2000\phantom{\rule{thinmathspace}{0ex}}{\mathrm{k}}_{\mathrm{B}}\mathrm{T}$, whereas large values of ${k}_{\mathrm{off}}$ result in $\mathrm{\Delta}{E}_{\mathrm{r}\mathrm{b}\mathrm{c}}\approx 100\phantom{\rule{thinmathspace}{0ex}}{\mathrm{k}}_{\mathrm{B}}\mathrm{T}$. The main reason is that low values of ${k}_{\mathrm{off}}$ lead to a significant increase in the lifetime of individual bonds, allowing the parasite to form more bonds and thereby increase its adhesion energy and induce larger membrane deformations. Similarly, large values of ${k}_{\mathrm{off}}$ decrease the bond lifetime, resulting in a decrease in the adhesion energy. For instance, in case of ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}=0.5$, the parasite forms on average about 200 bonds, whereas for ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}=4$, the average number of bonds is approximately 15 (see Figure 6—figure supplement 1). Furthermore, a larger onrate for the short bonds yields a slight increase in the strength of membrane deformations in comparison to a smaller ${k}_{\mathrm{on}}^{\mathrm{short}}$.
Figure 6b,c shows that there is a clear correlation between the level of membrane deformations and average alignment time. For example, for offrates ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}\le 2$, the alignment times are comparable with those for the reference parameter case, while for offrates ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}>2$, there is a strong increase in alignment times, which is correlated with insignificant membrane deformations. A shorter alignment time for ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}\le 2$ is due to the partial wrapping of the parasite by the RBC membrane, which is consistent with the previous study by Hillringhaus et al., 2019 that demonstrates the importance of membrane deformation for merozoite alignment. Note that the fixedtime displacement $\mathrm{\Delta}d$ in Figure 6a significantly increases with ${k}_{\mathrm{off}}$ due to a weaker adhesion. This seems to imply that the parasite alignment may proceed faster for ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}>2$. However, as it is evident from Figure 6c, this simple expectation is not applicable here, indicating that a faster motion of the parasite at the RBC surface may not necessarily result in a faster alignment. Alignment times for ${k}_{\mathrm{on}}^{\mathrm{short}}/{k}_{\mathrm{on}}^{\mathrm{long}}=8$ are generally shorter than for ${k}_{\mathrm{on}}^{\mathrm{short}}/{k}_{\mathrm{on}}^{\mathrm{long}}=4$ because of a slightly stronger parasite wrapping by the membrane. A seemingly opposite result for ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}=0.5$ in Figure 6c is likely due to insufficient statistics in the probability maps used for MC sampling, as they are constructed based on several direct simulations. Accurate resolution of small differences in alignment times is challenging, as it requires a large number of direct simulations.
Another bond parameter, which may affect parasite alignment, is the extensional rigidities of both bond types. Figure 7 presents RBC deformation energy and the number of bonds for five times softer and stiffer bonds than those in the reference case. Bonds with a larger rigidity lead to the formation of a larger number of bonds, more membrane wrapping, and a larger RBC deformation energy in comparison to soft bonds. The physical mechanism is that stiffer bonds facilitate a smaller distance between the membrane and the parasite at the edge of adhesion area between them, which favors further wrapping by the formation of additional bonds. Therefore, the spring rigidity in our model can mediate distancelimited bond formation at the edge of adhesion area between the parasite and the membrane, which affects merozoite alignment (see Figure 7—figure supplement 1), and is connected to membrane bending rigidity and the degree of wrapping. Consistently, simulations of the merozoite on a rigid RBC show no effect of the bond extensional rigidities on parasite alignment, because no significant membrane deformations are induced by parasite adhesion.
Furthermore, we consider effect of the density of long ligands ${\rho}_{\mathrm{long}}$ on parasite alignment. For the reference parameter set, ${\rho}_{\mathrm{long}}$ is chosen to be ${\rho}_{\mathrm{long}}/{\rho}_{\mathrm{para}}=0.4$, so that ${\rho}_{\mathrm{short}}/{\rho}_{\mathrm{para}}=0.6$. Figure 8 presents the number of short and long bonds as well as parasite alignment times as a function of ${\rho}_{\mathrm{long}}/{\rho}_{\mathrm{para}}$. Interestingly, the number of short bonds increases with increasing ${\rho}_{\mathrm{long}}$, even though the density of short ligands ${\rho}_{\mathrm{short}}$ decreases. This occurs due to the fact that more long bonds further stabilize parasite adhesion, allowing the formation of more short bonds. Note that for the density ${\rho}_{\mathrm{long}}/{\rho}_{\mathrm{para}}=0.1$ in Figure 8b, the value of $\u27e8{t}_{n}\u27e9$ is omitted, as the alignment criteria in Equation 4 have not successfully been met during the entire course of direct simulations, yielding the probability of parasite alignment in MC sampling to be zero. For ligand densities ${\rho}_{\mathrm{long}}/{\rho}_{\mathrm{para}}\ge 0.3$, both bond numbers and alignment times remain nearly independent of ${\rho}_{\mathrm{long}}$. However, the average alignment time for ${\rho}_{\mathrm{long}}/{\rho}_{\mathrm{para}}=0.2$ is about $30\mathrm{s}$ which is roughly three times longer than for the reference case. Note that $30\mathrm{s}$ is longer than the total length ($\approx 26\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}$) of direct simulations. Nevertheless, parasite alignment has occurred in some of these simulations, resulting in a small nonzero probability of merozoite alignment and a relatively long $\u27e8{t}_{n}\u27e9$ calculated through the MC sampling. The fact that $\u27e8{t}_{n}\u27e9$ for ${\rho}_{\mathrm{long}}/{\rho}_{\mathrm{para}}=0.2$ is longer than the total time of direct simulations means that the probability of parasite alignment is likely overestimated, indicating that the average alignment time should be even longer than $30\mathrm{s}$. An increase of $\u27e8{t}_{n}\u27e9$ with decreasing values of ${\rho}_{\mathrm{long}}$ is consistent with a significant decrease in membrane deformations (see Figure 8—figure supplement 1). For offrates ${k}_{\mathrm{off}}<72.6{\tau}^{1}$, the trends illustrated in Figure 8 remain qualitatively the same.
The importance of different ligand densities discussed above triggers the question whether both ligand types are necessary. Simulations performed with only short ligands (i.e., ${\rho}_{\mathrm{short}}/{\rho}_{\mathrm{para}}=1$) for several different ${k}_{\mathrm{off}}$ rates show that the parasite is not able to achieve significant wrapping by the membrane, because such ligands are too short to facilitate progressive membrane attachment over a curved parasite surface. This limitation is directly connected to the density of available receptors on the RBC surface, which is determined in our model by the membrane resolution. For the same reason, parasite mobility is impaired, as it is largely mediated by bond formation/dissociation at the edge of adhesion area between the parasite and the membrane. Therefore, the model with only short ligands does not reproduce proper parasite alignment. Simulations performed with only long ligands (i.e., ${\rho}_{\mathrm{long}}/{\rho}_{\mathrm{para}}=1$) show that the parasite mobility and alignment can be well reproduced, see Figure 8—figure supplement 2. Thus, the presence of long bonds aids in the stabilization of merozoite adhesion and the enhancement of parasite motion, such that long bonds serve as some sort of effective leverages. Theoretically, a model with only long ligands would be sufficient to reproduce the proper parasite alignment; however, current biomolecular knowledge about parasite coating does not support the presence of many bonds with a length of about $100\mathrm{nm}$. We speculate that short bonds are necessary (i) to stabilize parasite adhesion, as the density of long ligands is likely low, and (ii) to bring the two cells in sufficiently close contact (about $10\mathrm{nm}$) to facilitate the formation of a tight junction required for invasion. Thus, the presence of both ligand types is likely necessary for a successful invasion.
Effect of RBC rigidity
To investigate the effect of RBC rigidity on the alignment of a merozoite, we consider a nearly rigid cell membrane by increasing both bending rigidity and Young’s modulus by two orders of magnitude in comparison to a healthy RBC. Such a rigid RBC shows no significant membrane deformations for the reference interaction parameters given in Table 2, see Video 3. Comparison of parasite fixedtime displacements and alignment times for flexible and rigid membranes is shown in Figure 9 for two different values of ${k}_{\mathrm{off}}$. Clearly, larger RBC rigidity leads to much longer parasite alignment times (see Figure 9b), emphasizing again the importance of membrane deformations for merozoite alignment. For offrates ${k}_{\mathrm{off}}/{k}_{\mathrm{on}}^{\mathrm{long}}<2$, parasite alignment at the surface of a rigid RBC is not achieved within the course of the simulation. As the offrate increases, alignment time at the rigid membrane becomes comparable with that for the flexible membrane, because large enough ${k}_{\mathrm{off}}$ values do not result in strong membrane deformations even for the flexible RBC. Thus, for large offrates, the parasite’s alignment solely relies on its rotational dynamics controlled by the bond kinetic rates.
Figure 9a presents a comparison of parasite fixedtime displacements at the flexible and rigid membranes. In both cases, parasite displacements increase with increasing ${k}_{\mathrm{off}}$, as expected. However, the displacement at the rigid membrane is larger than at the flexible membrane (for visual comparison, see Videos 1 and 3), because the merozoite forms less bonds at the rigid surface. For the same reason, the variance of parasite displacements is larger for the rigid RBC than for the flexible RBC. Note that an increase in ${k}_{\mathrm{off}}$ results in an increase of fixedtime displacement and a decrease of alignment time for the rigid membrane, whereas for flexible RBC, an increase in offrate leads to an elevation of both fixedtime displacement and alignment time. This implies that for a rigid RBC, fast kinetics or weak adhesion are favorable for a quick alignment. In contrast, for a flexible RBC, slow kinetics or strong adhesion are advantageous for fast alignment, since the parasite employs RBC deformation for efficient alignment by partial membrane wrapping.
Discussion and conclusions
We have investigated the alignment of a merozoite at RBC membrane using a realistic twostate bonddynamics model for parasite adhesion. Motivated by experiments (Bannister et al., 1986), parasite adhesion is modeled by two bond types, with long and short binding ranges. Since RBCparasite interactions and the corresponding bond properties are experimentally not yet well characterized, the calibration of bond parameters is based on parasite fixedtime displacement at the membrane from existing experiments (Weiss et al., 2015), which is in the range of $0.30.8\phantom{\rule{thinmathspace}{0ex}}\mu m$. The presented model is able to reproduce quantitatively experimentally measured alignment times. Simulated alignment times are in the range between a few seconds and $26\mathrm{s}$, while the analysis of experimental videos by Weiss et al., 2015 yields an average alignment time of $16\mathrm{s}$. Another independent experimental study by Yahata et al., 2012 reports alignment times in the range between 7 and $44\mathrm{s}$, which agree relatively well with our simulation predictions. In addition to the good agreement between simulated and experimental alignment times, our model reproduces well dynamic RBC membrane deformations frequently observed in experiments (Dvorak et al., 1975; Gilson and Crabb, 2009; Crick et al., 2013).
Our main result is that parasite alignment is mediated by RBC membrane deformations and a diffusivelike dynamics due to the stochastic nature of parasitemembrane interactions. Average number of bonds $\u27e8{n}_{b}\u27e9$ between the parasite and the membrane is governed by the ratio ${k}_{\mathrm{on}}/{k}_{\mathrm{off}}=\mathrm{exp}(\mathrm{\Delta}{U}_{b}/{k}_{\mathrm{B}}T)$ that is connected to the binding energy $\mathrm{\Delta}{U}_{b}$ of a single bond and determines the strength of membrane deformations. Our results show that membrane deformations speed up the alignment through partial wrapping of the parasite, facilitating a contact between the parasite apex and the membrane. This conclusion is consistent with the previous simulation study (Hillringhaus et al., 2019), where merozoite adhesion has been modeled by a laterally homogeneous interaction potential whose strength controls RBC deformations. The importance of membrane deformation is also corroborated by simulations of parasite alignment at a rigid RBC, which show a drastic increase in alignment times. For a rigid membrane, the parasite alignment depends mainly on bond lifetime (i.e., ${\tau}_{b}\simeq 1/{k}_{\mathrm{off}}$), indicating that a low ${k}_{\mathrm{off}}$ or large bond lifetime may significantly decelerate the parasite’s rotational motion, and hence, increase its alignment time drastically. This conclusion agrees well with a recent simulation study (Jana and Mognetti, 2019) on the dynamics of two adhered colloids, whose effective rotational diffusion is governed not only by $\u27e8{n}_{b}\u27e9$ but also by ${\tau}_{b}$. Clearly, ${\tau}_{b}$ is also important for parasite dynamics at a deformable RBC, in addition to the membrane relaxation time ${\tau}_{p}$ on the scale of parasite size. The poor alignment of the merozoite at a stiff membrane can be a contributing factor, limiting parasite invasion. For example, infected RBCs in malaria become significantly stiffer than healthy cells (Suresh et al., 2005; Fedosov et al., 2011), limiting secondary invasion events. Furthermore, an increased RBC membrane stiffness is relevant in many other diseases, such as sickle cell anemia (Barabino et al., 2010), thalassemia (Peters et al., 2011), and stomatocytosis (Caulier et al., 2018), whose carriers are generally less susceptible to malaria infection.
For large values of ${k}_{\mathrm{off}}$, the parasite is not able to induce strong deformations even at a flexible membrane, so that the alignment times at rigid and deformable RBCs become comparable, and the alignment is governed solely by a diffusivelike rotational dynamics. The diffusivelike motion of the parasite at the membrane surface is facilitated by stochastic formation/dissociation of bonds between the two cell surfaces, and leads occasionally to a successful alignment. Therefore, our model is also able to explain the possibility of RBC invasion by a merozoite without preceding membrane deformations, which is observed much less frequently than the invasion preceded by significant RBC deformations (Weiss et al., 2015). Note that the RBCparasite adhesion model based on a laterally homogeneous interaction potential (Hillringhaus et al., 2019) predicts the complete failure of parasite alignment without significant membrane deformations, because it does not capture a diffusivelike rotational dynamics of the parasite. Thus, the bondbased model is more appropriate for the representation of RBCparasite interactions.
Even though the bond parameters in Table 2 were calibrated by the parasite fixedtime displacement obtained from experiments (Weiss et al., 2015), such a choice is likely not unique as some other set of parameters (e.g., receptor and ligand densities, bond rigidities and kinetic rates) may lead to statistically similar displacement characteristics. Nevertheless, it is important to emphasize that the discrete bonds in simulations should be thought of as ‘effective’ bonds, which likely represent a small cluster of real molecular bindings. Furthermore, since the parasite displacement is mainly controlled by the bond kinetics, this calibration procedure is rather robust in identifying an appropriate range of bond properties. Another important aspect of this model is the necessity of sufficiently long ligands and bonds to facilitate dynamic motion of the parasite at RBC surface. Simulations with only short ligands show that the parasite fails to induce significant wrapping by the membrane, leading to very little alignment success. Therefore, the long bonds serve as leverages for stable parasite adhesion and its motion at the membrane. Even though simulations with only long ligands indicate that a proper alignment can be achieved in this case, the existence of a dense population of long bonds has currently no support experimentally. Furthermore, we hypothesize that short enough bonds are necessary to enable the formation of a tight junction for parasite invasion, which requires a contact distance of about $10\mathrm{nm}$ between the two cells. Thus, our simulations suggest that both ligand types are likely necessary.
Electron microscopy images of adhered parasites (Bannister et al., 1986) suggest that the density of long bonds can be as low as 5  10%. However, the density of long ligands and bonds in our simulations is limited by the resolution of both the RBC and parasite to be larger than about 20%. A much finer membrane model would alleviate this limitation, but it would be prohibitively expensive computationally. Note that such heterogeneous receptorligand interactions exist in other biological systems as well. For example, during leukocyte binding in the microvasculature, both selectin and integrin molecules participate in adhesion and work synergistically, even though they have distinct functions (Ley et al., 2007). Furthermore, infected RBCs in malaria adhere to endothelial cells via two distinct receptors, ICAM1 and CD36, where binding with ICAM1 exhibits a catchlike bond, while the interaction with CD36 is a sliplike bond (Lim et al., 2017).
Several studies (Cowman et al., 2012; Dasgupta et al., 2014; Singh et al., 2010) about RBCparasite interactions hypothesize the existence of an adhesion gradient along the parasite body, which is expected to facilitate alignment. Based on the RBCparasite adhesion model with a laterally homogeneous interaction potential (Hillringhaus et al., 2019), it was shown that an adhesion gradient, where the potential strength increases toward the apex of a merozoite, generally accelerates parasite alignment. No definite conclusions about possible gradients can be made in the context of that model, because even in the case of no adhesion gradients, it predicts very short alignment times of about two orders of magnitude smaller than measured experimentally. An introduction of adhesion gradients in our bondbased interaction model leads qualitatively to the following conclusions: (i) Weak adhesion gradients do not significantly disturb the irregular motion of a parasite at RBC membrane, and have a negligible effect on the alignment. (ii) Strong adhesion gradients often result in a controlled direct reorientation of the parasite toward its apex, suppressing the irregular motion observed experimentally. These preliminary results do not permit a definite conclusion about the possible existence of adhesion gradients, as moderate adhesion gradients may exist and aid partially in the alignment process. Nevertheless, our model shows that adhesion gradients are not necessary, since the main parasite properties, such as dynamic motion and realistic alignment times, can be reproduced well by the bondbased model without adhesion gradients.
In conclusion, our model suggests that the parasite alignment can be explained by the passive compliance hypothesis (Introini et al., 2018; Hillringhaus et al., 2019), such that no additional active mechanisms or processes are necessary. Of course, this does not eliminate the possible existence of some active mechanisms, which may participate in the alignment process. Another limitation of many studies is that the parasite alignment is investigated under static (no flow) conditions, whereas in vivo, parasite alignment and invasion occur under a variety of blood flow conditions, including different flow stresses and flowinduced RBC deformations (Lanotte et al., 2016). Further experiments are needed to investigate RBCparasite interactions for realistic bloodflow scenarios. The bondbased model proposed here is expected to be useful for the quantification of such experimental studies and for a better understanding of RBCparasite adhesion under blood flow conditions.
Model and methods
Red blood cell model
The total potential energy of the RBC model is given by Fedosov et al., 2010a; Fedosov et al., 2010b
Here, the term ${U}_{\mathrm{sp}}$ represents the elasticity of spectrin network, which is attached to the back side of the lipid membrane. ${U}_{\mathrm{bend}}$ models the resistance of the lipid bilayer to bending. ${U}_{\mathrm{area}}$ and ${U}_{\mathrm{vol}}$ constrain the area and volume of RBC membrane, mimicking incompressibility of the lipid bilayer and the cytosol, respectively.
The elastic energy term ${U}_{\mathrm{sp}}$ is given by
where the first term is the attractive wormlike chain potential, while the second term corresponds to a repulsive potential with a strength ${\lambda}_{i}$. Furthermore, ${\mathrm{\ell}}_{i}$ is the length of the ith spring, ${p}_{i}$ is the persistence length, ${\mathrm{\ell}}_{i}^{\mathrm{max}}$ is the maximum extension, and ${x}_{i}={\mathrm{\ell}}_{i}/{\mathrm{\ell}}_{i}^{\mathrm{max}}$. The stressfree state of the elastic network is considered to be a biconcave RBC shape, such that initial lengths in the triangulation of this shape define equilibrium spring lengths ${l}_{i}^{0}$. For a regular hexagonal network, its twodimensional (2D) shear modulus µ can be derived in terms of model parameters as (Fedosov et al., 2010a; Fedosov et al., 2010b)
where $\overline{x}={\mathrm{\ell}}_{i}^{0}/{\mathrm{\ell}}_{i}^{\mathrm{max}}$ is a constant for all $i$. Thus, for given values of µ, $\overline{x}$, and ${\mathrm{\ell}}_{i}^{0}$, individual spring parameters ${p}_{i}$ and ${\lambda}_{i}$ are calculated by using Equation 8 and the force balance $\mathrm{\partial}{U}_{\mathrm{s}\mathrm{p}}/\mathrm{\partial}{l}_{i}{}_{{l}_{i}^{0}}=0$ for each spring.
The bending energy of the membrane is expressed as (Gompper and Kroll, 1996; Gompper and Kroll, 2004)
where $\kappa $ is the bending modulus, ${\mathbf{\mathbf{n}}}_{i}^{\mathrm{rbc}}$ is a unit normal of the membrane at vertex $i$, ${\sigma}_{i}=\left({\sum}_{j(i)}{\sigma}_{ij}{r}_{ij}\right)/4$ is the area of dual cell of vertex $i$, and ${\sigma}_{ij}={r}_{ij}[\mathrm{cot}({\theta}_{1})+\mathrm{cot}({\theta}_{2})]/2$ is the length of the bond in dual lattice, with the two angles ${\theta}_{1}$ and ${\theta}_{2}$ opposite to the shared bond ${\mathbf{\mathbf{r}}}_{ij}$.
The last two terms in Equation 6,
constrain surface area and volume of the RBC (Fedosov et al., 2010a; Fedosov et al., 2010b), where ${k}_{\mathrm{a}}$ and ${k}_{\mathrm{\ell}}$ control the total surface area $A$ and local areas ${A}_{i}$ of each triangle to be close to desired total area ${A}_{0}$ and local areas ${A}_{i}^{0}$, respectively. The coefficient ${k}_{\mathrm{v}}$ controls the total volume $V$ of the cell. The values of these coefficients are chosen large enough such that the area and volume fluctuate within 1% of the desired values.
The elasticity of a healthy RBC is characterized by the 2D shear modulus $\mu \approx 4.8\mu \mathrm{N}{\mathrm{m}}^{1}$, which corresponds to the 2D Young’s modulus $\mathrm{Y}\approx 18.9\mu \mathrm{N}{\mathrm{m}}^{1}$ for a nearly incompressible membrane (Suresh et al., 2005; Fedosov et al., 2010a). These values are employed in all simulations unless stated otherwise. The described membrane model has been shown to accurately capture RBC mechanics (Fedosov et al., 2010a; Fedosov et al., 2010b) and membrane fluctuations (Turlier et al., 2016).
RBCparasite adhesion interaction
Interaction between parasite and RBC membrane has two components. The first part imposes excludedvolume interactions between the RBC and merozoite (i.e. no overlap between them), using the purely repulsive part of the LennardJones (LJ) potential
This potential acts between every pair of RBC and parasite vertices separated by a distance $r=\left{\mathbf{\mathbf{r}}}_{\mathrm{rbc}}{\mathbf{\mathbf{r}}}_{\mathrm{para}}\right$ that is smaller than ${2}^{1/6}\sigma$. Here, $\u03f5$ represents the strength of interaction and $\sigma $ is the characteristic length scale of repulsion. The distance $\sigma $ can be thought of as an effective membrane thickness (imagine a surface constructed from overlapping spheres with a diameter $\sigma $). Normally, $\sigma $ should be selected as small as possible for a given resolution length of both the RBC membrane and parasite, which is about 0.2 μm in our models. Therefore, $\sigma =0.2\phantom{\rule{thinmathspace}{0ex}}\mu \text{m}$ is chosen, such that no overlap between the cells is guaranteed and the interacting surface is smooth enough.
The attractive part of RBCparasite interaction is modeled by a reversible twostate bond model. Bonds can form between RBC membrane vertices representing receptors and merozoite vertices corresponding to ligands, while existing bonds can also dissociate. These bonds represent RBCparasite adhesion through existing agonists at the surface of these cells and can be formed by two different types of ligands:
long ligands with an effective binding range ${\ell}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{long}=100\phantom{\rule{thinmathspace}{0ex}}nm$,
short ligands with an effective binding range ${\ell}_{eff}^{short}=20\phantom{\rule{thinmathspace}{0ex}}nm$,
which is motivated by electron microscopy observations of RBCmerozoite adhesion (Bannister et al., 1986). Long ligands result in long bonds, while short ligands lead to short bonds. Both bond types are modeled by harmonic springs with the potential energy given by
where ${\lambda}_{\mathrm{type}}$ is the spring extensional rigidity of either long or short bond type and ${\mathrm{\ell}}_{0}$ is the equilibrium bond length. To model the dynamic twostate interaction, constant (i.e. length independent) on and offrates (${k}_{\mathrm{on}}^{\mathrm{short}}$, ${k}_{\mathrm{on}}^{\mathrm{long}}$, and ${k}_{\mathrm{off}}$) are chosen, in order to simplify the model and reduce the number of parameters. Furthermore, the offrate for both bond types is selected to be same. Note that this model can easily be extended to lengthdependent rates.
To implement the different bond types, each vertex at the parasite surface represents either a long or a short ligand. The choice of vertices that correspond to long or short ligands is made randomly for fixed ligand densities ${\rho}_{\mathrm{long}}$ and ${\rho}_{\mathrm{short}}$. To avoid possible artifacts of a single discrete ligand distribution, each independent simulation assumes a different random choice of ligands with their respective densities kept constant. Bonds between the vertices at the RBC and parasite surfaces can form if the distance between two vertices is smaller than the corresponding cutoff distances ${\mathrm{\ell}}_{0}+{\mathrm{\ell}}_{\mathrm{eff}}^{\mathrm{long}}$ and ${\mathrm{\ell}}_{0}+{\mathrm{\ell}}_{\mathrm{eff}}^{\mathrm{short}}$, which remain the same in all simulations. Here, ${\ell}_{0}={2}^{1/6}\sigma$ corresponds to the length of the excludedvolume LJ interactions between the vertices of RBC and parasite, whose choice is defined by a characteristic discretization length of the RBC membrane. Only a single bond is allowed at each vertex for the both ligand types. Note that existing bonds can stretch beyond their effective binding ranges ${\mathrm{\ell}}_{\mathrm{eff}}^{\mathrm{long}}$ and ${\mathrm{\ell}}_{\mathrm{eff}}^{\mathrm{short}}$.
Hydrodynamic interactions
Hydrodynamic interactions are modeled using the dissipative particle dynamics (DPD) method (Hoogerbrugge and Koelman, 1992; Español and Warren, 1995), where fluid is represented by a collection of particles interacting through three types of pairwise forces: conservative ${\mathbf{\mathbf{F}}}_{ij}^{C}$, dissipative ${\mathbf{\mathbf{F}}}_{ij}^{D}$, and random ${\mathbf{\mathbf{F}}}_{ij}^{R}$ forces. The total force between particles $i$ and $j$ is given by
The conservative force models fluid compressibility, whereas the dissipative and random forces maintain a desired temperature of the system. The dissipative force also gives rise to fluid viscosity, which is generally measured in DPD by simulating a reversiblePoiseuille flow (Backer et al., 2005; Fedosov et al., 2010c). The DPD interactions are implemented only between the pairs of fluidfluid, fluidRBC, and fluidparasite particles. DPD interaction parameters are selected such that they impose noslip boundary condition at RBC and parasite surfaces (Fedosov et al., 2010a; Hillringhaus et al., 2019).
Simulation setup
Simulation domain with a size of $7.7{D}_{0}\times 3.1{D}_{0}\times 3.1{D}_{0}$ contains both RBC and parasite suspended in a DPD fluid, where ${D}_{0}=\sqrt{{A}_{0}/\pi}$ is the effective RBC diameter. Periodic boundary conditions are imposed in all directions. Initially, the parasite is placed close enough to the RBC membrane, so that the interaction between them is immediately possible. The initial parasite orientation is with its apex directed away from the membrane to mimic least favorable attachment configuration.
The main simulation parameters are shown in Table 1, both in simulation and physical units. To compare simulation units to physical units, a basic length scale is defined as the effective RBC diameter ${D}_{0}$, an energy scale as ${k}_{\mathrm{B}}T$, and a time scale as RBC membrane relaxation time $\tau =\eta {D}_{0}^{3}/\kappa $. For average properties of a healthy RBC, the effective diameter is $\mathrm{D}}_{0}\simeq 6.5\mu \mathrm{m$ with $\mathrm{A}}_{0}=133.5\mu {\mathrm{m}}^{2$ (Evans and Skalak, 1980) and the relaxation time becomes $\tau \approx 0.92\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}$ for the bending modulus $\kappa =3\times {10}^{19}\phantom{\rule{thinmathspace}{0ex}}J$ (Evans, 1983; Fedosov et al., 2010a) and plasma viscosity $\eta =1\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{P}\mathrm{a}\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}$. All simulations are performed on the supercomputer JURECA Jülich Supercomputing Centre, 2018 at the Jülich Supercomputing Centre, Forschungszentrum Jülich.
MonteCarlo sampling of alignment times
One of the main foci of our study is to obtain distributions of parasite alignment times for various conditions, which requires a large number of simulations of merozoite alignment. In order to significantly reduce the computational effort, MonteCarlo (MC) sampling of alignment times, which is guided by direct DPD simulations of RBCparasite adhesion, is employed. The MC sampling is based on a twodimensional probability map (see e.g. Figure 4a), which characterizes parasite orientation at the membrane surface through the distance $d}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{x}$ between the parasite apex and membrane and merozoite alignment angle $\theta$ (see Figure 3a for definitions of $d}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{x}$ and $\theta$). To construct such a probability map, possible $d}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{x}$ and $\theta$ values are binned into a number of orientation states $(i,j)=({d}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{x}}^{i},{\theta}^{j})$, and the probability $P(i,j)$ of each state is computed from at least 10 long DPD simulations of RBCparasite adhesion. We have verified that 10 independent DPD simulations are enough to reliably compute a probability map through its convergence with the number of DPD simulations. In the MC algorithm, changes in parasite orientation are modeled by transitions between different states, using the Metropolis algorithm. Thus, the transition from a state $(i,j)$ to one of the neighboring states $(i+1,j)$, $(i1,j)$, $(i,j+1)$ or $(i,j1)$ is selected randomly with a probability of $1/4$, and this move is accepted if $\xi <P(\mathrm{n}\mathrm{e}\mathrm{w}\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e})/P(i,j)$, where $\xi$ is a random number drawn from a uniform distribution in the interval $[0,1]$. In summary, the MC sampling algorithm is performed as follows:
Initial parasite orientation is selected randomly by choosing a state $({d}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{x}}^{i},{\theta}^{j})$, which has a nonzero probability.
Transitions between the neighboring states are modeled according to the Metropolis algorithm described above.
MC procedure is stopped whenever predefined alignment criteria are reached, and the number of MC steps is interpreted as alignment time.
Note that the MC sampling algorithm fulfills detailed balance, but does not account for hydrodynamic interactions. The fulfillment of detailed balance for the Metropolis algorithm in equilibrium means that changes between different states $(i,j)$ and $({i}^{\mathrm{\prime}},{j}^{\mathrm{\prime}})$ (with energies $E}_{(i,j)$ and $E}_{({i}^{\mathrm{\prime}},{j}^{\mathrm{\prime}})$, respectively) are performed according to transition rates proportional to $\mathrm{exp}\left[({E}_{(i,j)}{E}_{({i}^{\mathrm{\prime}},{j}^{\mathrm{\prime}})})/({k}_{B}T)\right]$, which are directly connected to probabilities of different states. Noteworthy, the MC sampling is a fast and efficient way to sample the distribution of parasite alignment times.
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files. Source data files have been provided for Figures 29, including all figure supplements. Figure 1 is a model schematic, which does not contain any data.
References

Poiseuille flow to measure the viscosity of particle model fluidsThe Journal of Chemical Physics 122:154503.https://doi.org/10.1063/1.1883163

Structure and development of the surface coat of erythrocytic merozoites of Plasmodium knowlesiCell and Tissue Research 245:281–290.https://doi.org/10.1007/BF00213933

Sickle cell biomechanicsAnnual Review of Biomedical Engineering 12:345–367.https://doi.org/10.1146/annurevbioeng070909105339

Merozoite surface proteins in red blood cell invasion, immunity and vaccines against malariaFEMS Microbiology Reviews 40:343–372.https://doi.org/10.1093/femsre/fuw001

Primary red cell hydration disorders: pathogenesis and diagnosisInternational Journal of Laboratory Hematology 40:68–73.https://doi.org/10.1111/ijlh.12820

The cellular and molecular basis for malaria parasite invasion of the human red blood cellThe Journal of Cell Biology 198:961–971.https://doi.org/10.1083/jcb.201206112

Statistical mechanics of dissipative particle dynamicsEurophysics Letters 30:191–196.https://doi.org/10.1209/02955075/30/4/001

A multiscale red blood cell model with accurate mechanics, rheology, and dynamicsBiophysical Journal 98:2215–2225.https://doi.org/10.1016/j.bpj.2010.02.002

Systematic coarsegraining of spectrinlevel red blood cell modelsComputer Methods in Applied Mechanics and Engineering 199:1937–1948.https://doi.org/10.1016/j.cma.2010.02.001

Steady shear rheometry of dissipative particle dynamics models of polymer fluids in reverse poiseuille flowThe Journal of Chemical Physics 132:144103.https://doi.org/10.1063/1.3366658

Multiscale modeling of blood flow: from single cells to blood rheologyBiomechanics and Modeling in Mechanobiology 13:239–258.https://doi.org/10.1007/s1023701304979

Morphology and kinetics of the three distinct phases of red blood cell invasion by Plasmodium falciparum merozoitesInternational Journal for Parasitology 39:91–96.https://doi.org/10.1016/j.ijpara.2008.09.007

Random surface discretizations and the renormalization of the bending rigidityJournal De Physique I 6:1305–1320.https://doi.org/10.1051/jp1:1996246

BookTriangulatedsurface models of fluctuating membranesIn: Nelson DR, Piran T, Weinberg S, editors. Statistical Mechanics of Membranes and Surfaces. Singapore: World Scientific. pp. 359–426.

BookRigid Body Mechanics: Mathematics, Physics and ApplicationsWeinheim: WileyVCH.

Importance of erythrocyte deformability for the alignment of malaria parasite upon invasionBiophysical Journal 117:1202–1214.https://doi.org/10.1016/j.bpj.2019.08.027

Anomalous transport in the crowded world of biological cellsReports on Progress in Physics 76:046602.https://doi.org/10.1088/00344885/76/4/046602

Simulating microscopic hydrodynamic phenomena with dissipative particle dynamicsEurophysics Letters 19:155–160.https://doi.org/10.1209/02955075/19/3/001

JURECA: modular supercomputer at Jülich Supercomputing CentreJounral of LargeScale Research Facilities 4:A132.https://doi.org/10.17815/jlsrf41211

Is invasion efficiency in malaria controlled by preinvasion events?Trends in Parasitology 23:481–484.https://doi.org/10.1016/j.pt.2007.08.001

Getting to the site of inflammation: the leukocyte adhesion cascade updatedNature Reviews Immunology 7:678–689.https://doi.org/10.1038/nri2156

The role of calcium in the invasion of human erythrocytes by Plasmodium falciparumMolecular and Biochemical Parasitology 50:317–323.https://doi.org/10.1016/01666851(92)90229D

A biological interpretation of transient anomalous subdiffusion. I. qualitative modelBiophysical Journal 92:1178–1191.https://doi.org/10.1529/biophysj.106.092619

Physics of adherent cellsReviews of Modern Physics 85:1327–1381.https://doi.org/10.1103/RevModPhys.85.1327
Article and author information
Author details
Funding
The authors declare that there was no funding for this work.
Acknowledgements
We would like to express our gratitude to Virgilio L Lew and Pietro Cicuta from the University of Cambridge for insightful discussions. Sebastian Hillringhaus acknowledges support by the International Helmholtz Research School of Biophysics and Soft Matter (IHRS BioSoft). We gratefully acknowledge the computing time granted through JARAHPC on the supercomputer JURECA (Jülich Supercomputing Centre, 2018) at Forschungszentrum Jülich.
Version history
 Received: February 29, 2020
 Accepted: May 17, 2020
 Accepted Manuscript published: May 18, 2020 (version 1)
 Version of Record published: June 3, 2020 (version 2)
Copyright
© 2020, Hillringhaus et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 723
 views

 137
 downloads

 9
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Physics of Living Systems
The paramount importance of mechanical forces in morphogenesis and embryogenesis is widely recognized, but understanding the mechanism at the cellular and molecular level remains challenging. Because of its simple internal organization, Caenorhabditis elegans is a rewarding system of study. As demonstrated experimentally, after an initial period of steady elongation driven by the actomyosin network, muscle contractions operate a quasiperiodic sequence of bending, rotation, and torsion, that leads to the final fourfold size of the embryos before hatching. How actomyosin and muscles contribute to embryonic elongation is investigated here theoretically. A filamentary elastic model that converts stimuli generated by biochemical signals in the tissue into driving forces, explains embryonic deformation under actin bundles and muscle activity, and dictates mechanisms of late elongation based on the effects of energy conversion and dissipation. We quantify this dynamic transformation by stretches applied to a cylindrical structure that mimics the body shape in finite elasticity, obtaining good agreement and understanding of both wildtype and mutant embryos at all stages.

 Physics of Living Systems
While inhomogeneous diffusivity has been identified as a ubiquitous feature of the cellular interior, its implications for particle mobility and concentration at different length scales remain largely unexplored. In this work, we use agentbased simulations of diffusion to investigate how heterogeneous diffusivity affects the movement and concentration of diffusing particles. We propose that a nonequilibrium mode of membraneless compartmentalization arising from the convergence of diffusive trajectories into lowdiffusive sinks, which we call ‘diffusive lensing,’ is relevant for living systems. Our work highlights the phenomenon of diffusive lensing as a potentially key driver of mesoscale dynamics in the cytoplasm, with possible farreaching implications for biochemical processes.