In silicoanalysis of the invasion mechanics and invasiveness of the plasmodium falciparum merozoite

Although there has been considerable progress in understanding the factors that determine the invasiveness of plasmodium falciparum merozoites, the collective role of the biophysical characteristics of erythrocyte deformability in the invasion process is poorly understood. Cell shape, cytoplasmic viscosity, and membrane stability are the main determinants of erythrocyte deformability, but it remains unknown how these properties affect the merozoite’s invasiveness. This study aimed to investigate computationally (i) the role of erythrocyte morphology and merozoite-induced erythrocyte membrane damage in merozoite invasion and (ii) the suitability of mechanical markers of merozoite-induced erythrocyte membrane damage for screening of invasion-blocking antimalarial drugs. Finite element models were developed to represent a human erythrocyte and a spherocyte, their invasion by a malaria merozoite, and erythrocyte compression and nanoindentation as mechanical assays for membrane damage. Smoothed particle hydrodynamics represented the erythrocyte cytoplasm, and merozoite-induced erythrocyte membrane damage was implemented with a constitutive model. The invasiveness of the merozoite decreases with increased erythrocyte sphericity associated with genetic disorders such as hereditary spherocytosis. The invasiveness is larger when membrane damage is induced in the erythrocyte at an early invasion stage than throughout the invasion process. The minimum force required for a malaria merozoite to invade a human erythrocyte was predicted to be 11 pN. The findings on the invasion mechanics can guide future studies into the invasiveness of the merozoite. The nanoindentation simulations point to the potential of nanoindentation to determine erythrocyte membrane damage for screening novel invasion-blocking anti-malaria drugs.

Erythrocyte; Hereditary spherocytosis; Malaria; Finite element method; Smoothed particle hydrodynamics; Drug screening (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 3 1. Introduction The invasion of malaria merozoites into human erythrocytes has been extensively studied as a potential target for antimalarial drugs (Flannery et al. 2013). During the invasion of erythrocytes, the merozoite is highly exposed to the host immune system and is highly vulnerable to therapeutic interventions; hence, it is an essential target for antimalarial drugs. Despite the enormous interest in studying the invasion process, fundamental gaps in the knowledge of the invasion process remain. To date, the impact of erythrocyte morphology and local erythrocyte membrane damage on the invasiveness of a merozoite has not been explored comprehensively due to the lack of detailed models for invasion mechanics to complement findings from experimental studies.

of in vivo approaches have limited the mechanistic understanding of the invasion process. In silico approaches provide alternative methods of ascertaining how biomechanical factors may contribute to the invasiveness of the merozoite. The merozoite entry into an erythrocyte is an active process that involves the application of actomyosin-based forces on the erythrocyte membrane. The forces are transmitted to the erythrocyte membrane through contact with the merozoite surface (Dasgupta et al. 2014). To date, a detailed analysis of the mechanistic role of the erythrocyte membrane and associated structure, i.e., the spectrin network involved in the invasion process, is limited to 2D analytical models (Abdalrahman and Franz 2017). Additionally, current analytical models do not incorporate the remodelling of the erythrocyte membrane and are limited to 10% of the invasion. Hence, there is a need to develop realistic 3D invasion models to account for all factors determining the merozoite’s invasiveness, i.e., merozoite-induced membrane damage. Hereditary spherocytosis is caused by genetic alteration of erythrocyte membrane proteins, leading to the formation of spherocytes. Previously, it has been documented that cells with (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 4 hereditary spherocytosis have abnormal protein structure and thus have a low susceptibility to infection by the merozoite (Eber and Lux 2004). Despite this finding, little is known about the higher invasion resistance of these cells. The current study aimed to computationally investigate the erythrocyte mechanics during malaria parasite invasion with emphasis on (i) the factors contributing to the merozoite’s invasiveness and (ii) the impact of local disruption of the spectrin network on the global mechanical properties of the erythrocyte for assessing the feasibility of mechanical markers for testing the efficacy of invasion blocking antimalarial drugs. 2. Materials and methods 2.1. Geometric modelling 2.1.1. Healthy erythrocyte The 3D biconcave geometry of the erythrocyte was described by: z = ± D0 ( 1 - 4 ( x2 + y2 ) D0 2 ) 1 2⁄ ( a0 + a1 x2 + y2 D0 2 + a2 ( x2 + y2 )2 D0 4 ) (1) with principal coordinate directions x, y, z, the diameter of the undeformed erythrocyte D0 = 7.82 µm, and shape parameters a0 = 0.0518, a1 = 2.026 and a2 = -4.491 (Figure 1 a). The generated erythrocyte model has a volume of 94.47 µm3 and a surface area of 135 µm2, consistent with the literature (Li et al. 2014). (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 5 Figure 1: a) Erythrocyte geometry, b) Spherocyte geometry, c) Dimension of Plasmodium Falciparum merozoite based on cryo-EM data from Dasgupta et al. (2014, Fig. 2), d) Geometry of rigid egg-shaped merozoite used in current study. 2.1.2. Spherocyte In individuals with spherocytosis, erythrocytes take a spherical form due to alterations of erythrocyte membrane proteins. The geometry was described with r2 = (x - kc)2 + (y - lc)2 + (z - nc)2 (2) where r = 2,6625 µm is the radius of the spherocyte (Li et al. 2016) and kc, lc, nc are centre coordinates of the spherocyte, where kc = lc = nc = 0. x, y, and z denote the coordinate points on the surface of the spherocyte (Figure 1 b). The surface area to volume ratio is 14% smaller for the spherical shape than the discoid geometry. 2.1.3. Plasmodium falciparum merozoite The merozoite shape has been previously described based on cryo-x-ray images of free merozoites (Dasgupta et al. 2014, Fig. 2). From these data, the mean physical dimensions of the merozoite were determined as follows: Length Lm = 1.98 ± 0.08 μm, width W = 1.40 ± 0.06 μm, volume Vactual = 1.71 ± 0.15 μm3 and surface area Aactual = 8.06 ± 0.72 μm2 (Figure 1 c). The 3D merozoite geometry (Figure 1 d) was generated with z (θ) = [2 Ra - Rb (1 - cos θ) ] (1 + cos θ )/4 (3) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 6 ρ (θ) = sin θ [2 Ra - Rb (1 - cos θ) ]/4 (4) where ρ and z are cartesian coordinates, Ra = 1 µm, Rb = 0.7 µm are parameters that determine the egg shape profile of the merozoite, and θ is a polar angle with 0° ≤ θ ≤ 180°. Images from electron microscopy, cryo-electron tomography, cryo-x-ray tomography and widefield deconvolution fluorescence imaging of a merozoite during invasion show negligible changes in merozoite shape throughout the invasion process (Zuccala et al. 2016). Hence, this study treated the merozoite geometry as a rigid body. 2.1.4. Tight junction between merozoite and erythrocyte membrane The merozoite pulls itself into the erythrocyte through the tight junction complexes, which it establishes after forming the invasion pit (Pinder et al. 2000; Preiser et al. 2000; Riglar et al. 2011). The tight junction was modelled as a deforming mechanical link between the merozoite surface and the erythrocyte membrane, forming an annulus-like structure to facilitate erythrocyte membrane wrapping. The annulus is defined as a circular ring with an internal diameter of 0.76 µm and a cross-sectional radius of 0.04 µm (Figure 1 d). 2.2. Constitutive modelling 2.2.1. Erythrocyte membrane with merozoite-induced damage During merozoite invasion, the erythrocyte membrane deformation was considered as mainly due to the mechanical loads exerted by the merozoite’s actomyosin machinery and other external sources, such as blood pressure, whereas the entropic deformation was considered negligible. Hence the Helmholtz free energy function for the erythrocyte membrane deformation was only represented as internal strain energy. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 7 The damage induced by the merozoite was modelled by modifying the strain energy density function of the erythrocyte membrane. Since the erythrocyte membrane comprises primarily an elastic spectrin network and can be considered an elastic, isotropic, and nearly incompressible continuum, the strain energy density function is usually presented in a decoupled form comprising deviatoric and isochoric terms (Li 2016). The combination of incompressibility and large deformation of a nearly incompressible hyperelastic material presents difficulties for a displacement-based finite element method as the constraint J = det F = 1 on the deformation field is highly nonlinear (Weiss 1994). To overcome this challenge, a displacement-based finite element scheme must invoke a small change measure of volumetric deformation. Consequentially, the deformation gradient was decomposed into the dilatational and deviatoric parts to apply separate numerical treatments to each part (Weiss 1994). Therefore, the deformation gradient F (Gilson and Crabb 2009)and the left Cauchy-Green strain tensor B were divided into the volume-changing (dilatational) and the volume- preserving (distortional) parts, an approach often used in elasto-plasticity (Ogden 1978). The strain energy density function of the isotropic erythrocyte membrane was expressed in terms of the left Cauchy deformation tensor as B = F · FT. (5) With F = R U and FT = RT UT, where R is a rotation matrix, and U is a stretch tensor, Eqn. (5) becomes B = RU · UT 𝐑T = U · UT, (6) showing that the left Cauchy deformation tensor B is a stretch tensor and isotropic. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 8 Hence the strain energy density function ψ of a damaged erythrocyte membrane can be written in terms of invariants of the left Cauchy-Green deformation tensor: ψ (I1, I2, J, τ ) = Rd [ ψ0(I1,I2) + ψ0 (J) ] (7) with I1 = tr(B), I2 = 1 2 (I1 2 - tr(B𝟐)), and J = √det (B), where B is the left Cauchy-Green tensor, Rd is the damage parameter, ψ0 is the strain energy density of an intact erythrocyte membrane, and τ is the indentation time in seconds. A continuum damage mechanics framework for material stiffness deterioration suitable for implementation in Abaqus Explicit was used to simulate chemical damage induced by a merozoite during the entry process. The model borrows concepts from various strain-based damage models for soft biological materials and biodegradable polymers where constitutive hydrolytic degradation and time-dependent behaviour are described (Vieira et al. 2011). Rd is a function of chemical and mechanical damage parameters defined as: Rd(B, J, τ) = [ 1 - β2] e-( β1 τ + β2 ψ0 ) (8) where B and J denote the left Cauchy green tensor and the volumetric strain of the intact erythrocyte membrane, respectively, τ is the simulation time, and ψ0 is the strain energy density per reference volume of the intact erythrocyte membrane, β1 is the chemical damage parameter, and β2 is the mechanical damage parameter. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 9 For β2 = 0, the damage mode is purely chemical, i.e., damage to the erythrocyte membrane is not due to deformation. Thus, β1 represents chemical damage due to various modes of phosphorylation, and β2 represents mechanical damage associated with the tearing or rupturing of the protein chains in the erythrocyte membrane skeleton. The erythrocyte's membrane near incompressibility was defined with a Poisson’s ratio ν = 0.499. For incompressible materials, the contribution of the volumetric strain energy density function is neglected since J = 1. Since the erythrocyte membrane is nearly incompressible, the deviatoric strain invariants are I̅1 = J-2/3 I1 (9) I̅2 = J-4/3 I2. (10) The variation of the strain energy potential is by definition equal to the internal virtual work per reference volume V0, and this can be written as: δWi= ∫ J ( S : δe - pδϵvol ) dV0 = ∫ δψ dV0 (11) where δψ = 2 [ ( ∂ψ ∂I̅1 + I̅1 ∂ψ ∂I̅2 ) B* - ∂ψ ∂I̅2 B* ∙ B* ] : δe + J ∂ψ ∂J δϵvol δe = δD - 1 3 δϵ vol I δD = sym ( δL ) = 1 2 ( δL + δLT ) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 10 δϵvol = I : δD δL = ∂δu ∂X S = σ + pI p = - 1 3 I : σ = - ∂ψ ∂J . Hence the deviatoric stress with damage can be rewritten as: S = 2Rd J DEV [ ( ∂ψ ∂I̅1 + I̅1 ∂ψ ∂I̅2 ) B* - ∂ψ ∂I̅2 B* ∙B* ] (12) S = 2Rd J [ ( ∂ψ ∂I̅1 + I̅1 ∂ψ ∂I̅2 ) (B* - trace ( B* ) 3 ) - ∂ψ ∂I̅2 ( B*∙ B* - trace ( B*∙ B* ) 3 ) ]. (13) Therefore, the stress (deviatoric stress and volumetric stress) with damage can be written as σij= Rd J [ 1 J2/3 ( ∂ ψ0 ∂ I̅1 + I̅1 ∂ ψ0 ∂ I̅2 ) Bij - (I̅1 ∂ ψ0 ∂ I̅1 +2I̅2 ∂ ψ0 ∂I̅2 ) δij 3 - 1 J 4 3 ∂ ψ0 ∂I̅2 BikBkj] + Rd∂ ψ0 ∂J δij (14) and further as σij= e-(β1τ) J [ 1 J2/3 (C10 + I̅1C01)Bij - (I̅1C10+2I̅2C01)δij 3 - 1 J 4 3 C01BikBkj+ K0inJ J δij], (15) where δ is the Kronecker delta function. ψ0 is the Mooney Rivlin strain energy density function per unit reference volume of the intact erythrocyte membrane. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 11 The developed erythrocyte membrane damage (EMD) model, Eqn. (15), to induce localised damage in the erythrocyte membrane was implemented using the VUMAT subroutine, whereas the Mooney Rivlin law model to describe the constitutive response of the intact erythrocyte membrane was implemented using the Abaqus materials module. The erythrocyte membrane damage model was verified using a single shell element model subjected to an equi-biaxial strain of 1.1. The verification involved comparing the true stress obtained with the VUMAT subroutine and the built-in Mooney Rivlin law for an intact erythrocyte membrane with the chemical damage parameter β1 = 0. For this case, the constitutive responses of the developed VUMAT subroutine and the built-in Mooney Rivlin law are expected to be identical. Thereafter, the single shell element model was used with various degrees of chemical damage to evaluate the stability of the developed erythrocyte membrane damage model using Drucker’s stability criterion. The material's relative compressibility also determines the erythrocyte membrane's mechanical response. The relative compressibility is the ratio of the initial bulk modulus K0 to the initial shear modulus µ0 of the material: Rc = K0 μ0 (16) where µ0 and K0 are defined by Eqns. (17) and (18), respectively, and large Rc values show that the material is less compressible: μ0 = 2 (C10 + C01) (17) K0 = 2 D1 (18) The Poisson’s ratio ν for hyperelastic materials is related to R by: (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 12 ν = 3Rc-2 6Rc+2 (19) Furthermore, D1 is expressed as: D1 = 3 (1 - 2 ν) 2 (C10 + C01) (1 + ν) (20) For this study, the Poisson’s ratio of the erythrocyte membrane was set to 0.499 to avoid numerical singularity; hence, D1 was set to 12 mm2/N. The material parameters for the Mooney Rivlin model are computed from the elastic modulus by: C10 = E 6 (1 + β) (21) and C01 = β C10 (Zhang and Zhang 2011), with β ranging from 0 to 0.5. For an elastic modulus E of 1 kPa and setting β = 0.1, the corresponding values of C10 and C01 are 152 Pa and 15.2 Pa, respectively. 2.2.2. Tight junction The tight junction was represented with an annulus structure, mimicking the function of the tight junction during the invasion process. The energy associated with the work done by the tight junction is not yet known. However, an estimate of the minimum energy contribution of the tight junction required for a successful invasion was determined with the developed finite element model, and the Mooney Rivlin law was used to define the mechanical response of the annulus structure. The material parameters were determined by minimising the objective function defined by: Fm (Cj) = ψm (Cj) × Dm (Cj) (22) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 13 where ψm is the strain energy density function for the annulus structure, and Dm is the diameter of the tight junction at maximum indentation depth at the end of the invasion when the simulation step time is 1.1 s. The binary search algorithm below was used to search for

parameters that minimise the objective function in Eqn. (22). Given an array Cj of n elements C0j, C1j, C2j …. Cnj-1 such that C0j ≤ C1j ≤ C2j ≤ … ≤ Cnj-1, the following pseudo-code uses the binary to find the value of the material parameter in Cj that minimises the objective function. 1. Set Lj = 0 and Rj = nj - 1 2. If Lj > Rj, the search terminates as unsuccessful 3. Determine the middle element index mj = floor ( [Lj + Rj ]/2 ) 4. Compute the Dm using the middle element parameter in Cj. 5. If Dm > Md + ē update Cj such that Cj ≥ Cm, compute Fm and go to step 2. 6. If Dm < Md + ē, update Cj such that Cj ≤ Cm, compute Fm and go to step 2. 7. If Dm ≈ Md the search is done, and compute Fmin = min (Fm) Here, Md is the width of the merozoite at maximum indentation depth, and ē is the clearance between the tight junction and the merozoite at the maximum indentation depth. When Fm = Fmin = min (Fm) (Figure 2 a), Cj gives a minimum strain energy ψm (Cj) such that Dm ≈ Md (Figure 2 b). The strain energy ψm (Cj) of the annulus structure increases from point a to b and c (Figure 2 a) while the diameter Dm decreases from point a to b but remains constant from point b to c (Figure 2 a and b). The annulus structure's mechanical properties mimicking the tight junction's role were determined by minimising the objective function, resulting in values of the Mooney Rivlin parameter of C10 = 0.04 MPa, C01 = 0.004 MPa and D1 = 0.3 mm2/N. The determined (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 14 mechanical properties represent the minimum energy the annulus structure requires to ensure erythrocyte membrane wrapping. Figure 2: Illustration of a) the minimisation process of the objective function and b) the determination of the Mooney Rivlin parameters. 2.2.3. Erythrocyte cytoplasm The smoothed particle hydrodynamics (SPH) method was used to model the erythrocyte cytoplasm deformation. Smoothed particle hydrodynamics is a fully Lagrangian mesh-free modelling scheme permitting the discretisation of a prescribed set of continuum equations by interpolating the properties directly at a discrete set of points distributed over the solution domain. This approach was first developed to solve PDE problems in astrophysics (Gingold and Monaghan 1977). In Abaqus, the SPH scheme discretises the continuum partial differential equations (Violeau and Rogers 2016). SPH uses an evolving interpolation scheme to approximate a field variable at any point in a domain. Using the particle approximation or (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 15 field function Eqn. (23) and its derivative Eqn. (24), the Navier-Stokes equation is discretised and solved using the explicit time integration method. fi = ∑ mj ρj N j=1 fj w ( | ri - rj | ,h ) = ∑ mj ρj N j=1 fj wij (23) ∇∙fi = ∑ mj ρj N j=1 fj ∙ ∇w ( | ri - rj |,h ) = ∑ mj ρj N j=1 fj ∙ ∇𝑤ij (24) where N is the total number of particles, h is the smoothing length, and ri and rj are the position vectors of the particle of interest and the particle in the neighbouring region, respectively. The field function f and its derivative are constructed using a smoothing or kernel function w (Wang et al. 2016, Fig. 1). Thus, the value of a variable at a particle of interest can be approximated by summing the contributions from a set of neighbouring particles, denoted by subscript j, for which the “kernel” function, w, is not zero. In Abaqus, the erythrocyte cytoplasmic domain is converted to SPH particles by activating a built-in conversion to SPH particle functionality in the mesh module. The erythrocyte cytoplasm primarily comprises viscous haemoglobin, mathematically described by the Navier-Stokes equation in the Lagrangian form (Ye et al. 2016). The erythrocyte cytoplasm is generally considered an incompressible Newtonian liquid, and thus, its dynamics are predicted by using the Navier-Stokes equations given by: ∇ ∙ v = 0 (25) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 16 dv dt = 1 ρ ( - ∇p + μ∇2 v ) +fext (26) where p, v, ρ, µ and fext represent the pressure, velocity, density, dynamic viscosity and external force vector, respectively. SPH solves the Navier-Stokes equations by discretising the whole computational domain into a set of particles. The Mie-Grüneisen equation of state, Eqn. (27), was used to model incompressible viscous laminar flow governed by the Navier- Stokes equation of motion. The volumetric response is governed by the equations of state, where the bulk modulus acts as a penalty parameter for the incompressible constraint. Since the viscosity of the erythrocyte cytoplasm is small, a small amount of shear resistance was specified in the materials module to suppress shear modes that can otherwise tangle the mesh. Here, the shear stiffness or viscosity was used as a penalty parameter. The default hourglass control was used because when the shear model is defined, the hourglass control forces are calculated based on the shear resistance of the erythrocyte cytoplasm, which provides very low shear strength, insufficient to prevent spurious hourglass modes. An equation of state is necessary for the erythrocyte cytoplasmic domain to link pressure P and density ρ. The Mie- Grüneisen equation of state used for this purpose (Monaghan 1988) is given as follows: p = ρ0c2ηm (1-ηms)2 (1- ηmΓ0 2 ) + ρ0Γ0Em (27) where ρ0 is the reference density, and c is the speed of sound, Γ0 = 0 is a material parameter, ηm = 1 − ρ0/ρ is the nominal volumetric compressive strain, Em is internal energy per unit mass. The background pressure P0 is added to avoid negative pressure values. The density is estimated from the particle distribution utilising the SPH interpolation. c and s define the linear relationship between the shock wave velocity, Us, and the particle velocity, Up, as follows: (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 17 Us = c + sUp (28) where s was set to zero such that Us = c =1,000 mm/s. 2.3. Finite element meshes The finite element model comprised multiple components with various element lengths, i.e., the erythrocyte membrane, the tight junction, the rigid merozoite, and the cytoplasm (Figure 3 a to d). The tight junction contains the smallest element length relative to the erythrocyte membrane and is generally used to determine stable time increments for the whole model. However, it was necessary to evaluate the stability of the solution due to the mesh density of the erythrocyte membrane. A mesh density study for the erythrocyte membrane was performed to determine the mesh density, which gives a stable solution. In Abaqus, the erythrocyte membrane was modelled using shell elements with a thickness of 0.01 µm (Hochmuth et al. 1973). SPH was applied for the erythrocyte cytoplasmic domain; hence, no mesh density study was performed for the cytoplasm. However, to ensure accurate simulation of the cytoplasm by using the SPH functionality, Abaqus requires more than 10,000 particles to be generated. Hence, 18,000 particles were generated using the particle conversion functionality. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 18 Figure 3: Finite element meshes and boundary conditions. Finite element mesh of a) erythrocyte membrane with triangular shell elements, b) erythrocyte cytoplasm with 8-node linear brick elements with reduced integration and hourglass control (C3D8R), c) annulus structure of the tight junction with ten-node modified quadratic tetrahedron elements (C3D10M), and d) rigid merozoite with three-node 3D rigid triangular elements (R3D3). Boundary conditions: e) Blood pressure applied on the outer surface of the erythrocyte membrane, f) the initial velocity of 2×10-30 mm/s applied at each node of the erythrocyte membrane, g) each node of the annulus structure is only allowed to displace in z and x directions, and h) the rigid merozoite is displaced in the negative y-direction. The erythrocyte membrane was meshed with 19,942 three-node triangular shell elements with reduced time integration (S3R). The reduced time integration algorithm provided more accurate results and reduced the running time. Conventional shell elements were preferred (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 19 because, unlike continuum shell elements, which only have displacement degrees of freedom, conventional shell elements have both displacement and rotational degrees of freedom. For the erythrocyte membrane to wrap around the merozoite effectively, each erythrocyte membrane node must have both displacement and rotation degrees of freedom. Erythrocytes undergo extreme deformations as they circulate through narrow capillaries in the human body. To accurately predict erythrocyte deformation, displacement and rotational degrees of freedom must be allowed in each element. Additionally, shell elements were preferred over 3D solid elements for the erythrocyte membrane because shell elements allow the modelling of thin features with fewer elements, thus reducing computational time. Shell elements are also easier to mesh and less prone to negative Jacobian errors, representing negative element volume or distorted elements that might occur when using extremely thin solid features. The erythrocyte cytoplasm was meshed using 841 eight-node linear brick elements with reduced integration and hourglass control (C3D8R). The mesh provides the initial spatial particle discretisation required for the SPH scheme. Seven hundred three-node 3D rigid triangular elements (R3D3) were used for the rigid merozoite, and 5,535 ten-node modified quadratic tetrahedron elements (C3D10M) were used for the annulus structure of the tight junction. 2.4. Boundary conditions Several assumptions and boundary conditions were considered: • The erythrocyte was assumed to be suspended in an Euclidean space with an initial velocity of approximately zero for each node . The load module was used to predefine the initial velocity fields at each node in the erythrocyte membrane (Figure 3 e). • A constant external pressure load of 16 kPa equal to systolic blood pressure was applied on the external surface of the erythrocyte membrane (Figure 3 f). (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 20 • The merozoite is displace d by 2 µm along its longitudinal axis of the merozoite. The displacement is applied using a ramp function from 0.1 s to 1.1 s (Figure 3 g). • Each node of the tight junction freely deforms in the x and z directions (Figure 3 h). 2.5. Contact interactions Two algorithms were used to model contact interactions between structures involved in the invasion process. The general contact algorithm was used to define contact interactions between the erythrocyte membrane, the erythrocyte cytoplasm and the tight junction. The contact pair algorithm was used to define contact between the merozoite and the outer surface of the erythrocyte membrane. The contact pair algorithm in Abaqus Explicit includes the contact surface weighting (balanced or pure master-slave) and the sliding formulation (finite, small, or infinitesimal). The contact pair algorithm with pure master-slave weighting was used for contact between the merozoite surface and the region of entry (ROE) on the erythrocyte membrane. The interacting surfaces can penetrate each other for the pure master- slave scheme, leading to numerical instabilities. To avoid numerical errors due to penetration, the mesh density of the slave surface must be greater than that of the master surface. A mesh density study for the two surfaces was conducted to determine the mesh sizes for both the master and slave surfaces. 2.6. Model validation 2.6.1. Validation of erythrocyte finite element model with simulation of a healthy erythrocyte in an optical tweezer/trap This section presents a detailed validation method of the developed erythrocyte model based on laser trap experimental data (Mills et al. 2004). The mechanical response of the erythrocyte predicted by the developed model is compared to data from an optical trap experiment during which a force of 193 pN is applied to stretch the erythrocyte (Mills et al. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 21 2004, Fig. 7). During an optical trap experiment, a stretching force is applied to one of the two silica beads attached to the erythrocyte. One microbead is fixed to the glass slide while the other is trapped using a laser beam (Song et al. 2017). Cell stretching is performed by moving the trapped microbead. The deformation is determined from images of undeformed and stretched erythrocytes. Similarly, during the finite element simulation of the optical trap experiment, the erythrocyte model was deformed in the axial and transverse directions due to a resultant force F = 200 pN applied on one microbead. 2.6.2. Validation of merozoite invasion finite element model Recently, Geoghegan et al. (2021) used a high spatiotemporal resolution lattice light-sheet microscopy (LLSM) to analyse the merozoite invasion into the erythrocyte by segmenting and tracking the formation of parasitophorous vacuole membrane (PVM). With this technique, the authors determined the portion of the erythrocyte membrane surface area that wraps the merozoite and the surface area of the erythrocyte membrane that does not wrap the merozoite. The study documents the decrease in the erythrocyte membrane surface area since a portion of the area wraps the merozoite, see Geoghegan et al. (2021, Fig. 2c). The data obtained from this study was used to validate the developed invasion model. In vivo, erythrocytes are subjected to physiological blood pressure. Hence, blood pressure must be applied to the surface of the erythrocyte model to obtain realistic simulation results with the invasion model. However, no physiological pressure was applied in the Geoghegan et al. (2021) experiment. As such, surface area data of the erythrocyte obtained from the invasion model where no surface pressure was applied was compared with erythrocyte areal data from the Geoghegan et al. (2021) experiment to validate the invasion model. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 22 2.7. Finite element analysis and case studies 2.7.1. Generalised explicit finite element analysis in Abaqus A generalised Abaqus Explicit dynamic analysis procedure was used to simulate the deformation of an erythrocyte during the invasion by a merozoite. This procedure involves numerically solving the momentum equilibrium Eqn. (29) using an explicit central difference time integration rule described by Eqns. (30) and (31). The momentum equilibrium is: MNJ ü i N = Pi J - Ii J (29) where MNJ is the mass matrix, PJ is the applied load vector, IJ is the internal force vector, uN denotes spatial degrees of freedom, and i is the increment number in the explicit dynamic step. MNJ is diagonalised to form a lumped mass matrix to reduce the computational complexity of the central difference time integration algorithm (Systemes 2015). u̇ i + 1/2 N = u̇ i - 1/2 N + ∆ti + 1 + ∆ti 2 ü i N (30) ui + 1 N = ui N + ∆ti + 1 u̇ i + 1/2 N (31) The dynamic explicit analysis was specified using the step module, where mass scaling and geometrically nonlinear analysis were used. Mass scaling was only applied to the erythrocyte membrane and the tight junction to obtain a quasi-static response. The quasi-static analysis was achieved by simulating the invasion process in the shortest time, i.e. for a period of 1.1 s, while keeping the inertial forces relatively low. The semi-automatic mass scaling was used throughout the step to scale mass elements periodically and effectively reduce the wave speed. The effectiveness of the mass scaling algorithm in ensuring a quasi-static solution was determined by ensuring that the total kinetic energy of the erythrocyte model is much smaller than the internal or strain energy of the erythrocyte model. Kinetic energy represents the effects of inertia on the global response of the erythrocyte model, while internal energy (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 23 represents the static effects. Time incrementation was achieved automatically by using built- in functionality. The adaptive, global estimation algorithm was applied to determine the maximum frequency of the entire model using the current dilatational wave speed. 2.7.2. Case Study 1: Impact of the erythrocyte morphology on the merozoite invasiveness With the developed finite element model, the impact of erythrocyte morphology and merozoite-induced damage on the invasiveness of the merozoite was assessed. The invasiveness of the merozoite in the convex region (Figure 4 a) and the concave region (Figure 4 b) was assessed by comparing the total invasion energies associated with the merozoite’s invasion of the convex and concave regions. To assess the impact of sphericity, i.e., surface area to volume ratio, the total invasion energies of two simulations are compared: One involving the entry of the merozoite into a normal discoid-shaped erythrocyte, and the other involving entry simulation of a merozoite into a spherical erythrocyte with the same volume as the normal erythrocyte (Figure 4 a). To date, it is unknown whether the spherical shape impacts the invasiveness of the merozoite. The morphometric parameters used to define the two erythrocyte morphologies, i.e., discoid and spherical morphologies, are given in Table 1. The developed spherical model of the erythrocyte represents a 27% reduction of a healthy erythrocyte’s surface-to-volume ratio S/V and total surface area. With the two models of the erythrocyte, i.e., spherical and discoid models, the impact of morphological variations on the invasiveness of the malaria parasite was investigated. Table 1: Erythrocyte parameters for case study 1. Erythrocyte shape Surface area S (µm2) Volume V (µm3) S/V (µm-1) Discoid 135 94 1.44 Spherical 98.47 94 1.05 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 24 Figure 4: a) Case study 1: Entry point of the merozoite (grey) into a human erythrocyte (blue) in a convex (left) and concave region (middle) and into a human spherocyte (blue) (right). b) Case study 2: Location and dimension of damage (red) induced in erythrocyte membrane. e) Case study 3: Initial configuration of the compression test simulation (left) and locations of induced membrane damage in concave and convex region of the erythrocyte. d) Case study 4: (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 25 Initial configuration of the nanoindentation simulation and locations of induced membrane damage in concave and convex region of the erythrocyte. 2.7.3. Case Study 2: Impact of phosphorylation induced damage in the erythrocyte membrane on the merozoite invasiveness Using the developed erythrocyte membrane damage model, two stages of erythrocyte membrane damage were investigated to assess the impact of the damage stage on the invasiveness of the merozoite. The first stage of damage was induced during the early invasion, i.e., from 0 to 0.1s, accounting for 9% of the invasion period. The second stage of damage involved damage induced throughout the entire invasion process. Erythrocyte membrane damage was induced in the red region, with a diameter of 1.028 µm (Figure 4 b). For each damage stage, the invasiveness of the merozoite was determined by evaluating the indentation forces induced by the merozoite. The damage parameters used were β1 = 5.4 and 11 for damage in the early invasion stage and β1 = 0.49 and 1 for damage throughout the invasion process. 2.7.4. Case Study 3: Erythrocyte compression The compression simulation with intact erythrocyte and erythrocyte with local membrane damage investigated the impact of local erythrocyte membrane damage on the erythrocyte's global mechanical and structural properties (Figure 4 c). First, an intact erythrocyte model was compressed to extract the compression data. Secondly, an erythrocyte with a locally damaged erythrocyte membrane was compressed to extract two compression data sets, one obtained by inducing damage in the concave region, i.e., damage location 1 and the other in the convex region of the erythrocyte membrane, i.e., damage location 2 (Figure 4 c). The diameter of the damage region was 1.028 µm. The compression plate (10 × 10 µm) was displaced by 1.3 µm for a period of 1.1 s, while the support plate (10 × 10 µm) was fixed to facilitate compression of the erythrocyte. The compression and the support plates were (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 26 modelled using 5,000 linear triangular rigid shell elements (R3D3). The contact pair algorithm was used to define contact between the compression plate-erythrocyte interface and the erythrocyte-support plate interface, where surface-to-surface contact formulation with kinematic contact method was defined. The normal and the tangential behaviour of the interfaces mentioned above were defined for this type of interaction. Isotropic tangential interaction with negligible friction coefficient was used to describe the tangential behaviour between the interfaces. Hard contact formulation, which allows separation of the interfaces mentioned above, was used to define the normal behaviour of the interfaces. The force- compression data obtained from indenting a locally damaged and intact erythrocyte were compared to assess the sensitivity of the global indentation for local erythrocyte membrane damage. Sufficient sensitivity may indicate the potential of compression tests for further development and implementation to study merozoite-induced damage in a physical experiment. The impact of surface pressure on the mechanical response of the locally damaged erythrocyte membrane was also investigated. 2.7.5. Case Study 4: Erythrocyte nanoindentation The nanoindentation simulation was conducted using a similar approach to the compression test simulation from the previous section. The nanoindentation simulation involved an intact erythrocyte and an erythrocyte with a locally damaged membrane to investigate the impact of local erythrocyte membrane damage (Figure 4 d). The spherical indenter with a diameter of 2 µm, representing 25.6% of erythrocyte diameter, was displaced by 1.31 µm, representing 64.7% of the maximum erythrocyte thickness. The support plate (10 × 10 µm) was fixed to facilitate indentation of the erythrocyte. The spherical indenter was modelled using 854 linear rigid triangular shell elements (R3D3), while the support plate was modelled by using 5,000 linear rigid triangular shell elements (R3D3). The contact pair algorithm was used to define contact between the spherical indenter-erythrocyte interface and the erythrocyte-support plate (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 27 interface, where the surface-to-surface contact formulation with kinematic contact method was defined. Both normal and tangential behaviour of the interfaces mentioned above were defined for this type of interaction. Isotropic tangential interaction with negligible friction coefficient was used to describe the tangential behaviour between the interfaces. Hard contact formulation that allows separation of the interfaces mentioned above was used to define the interfaces' normal behaviour. The force-indentation data obtained from the locally damaged and intact erythrocyte were compared to assess the sensitivity of the nanoindentation to local erythrocyte membrane damage. A sizeable difference between the force-indentation curves of the undamaged and damaged erythrocyte suggests that the nanoindentation test is sensitive to erythrocyte membrane damage. 3. Results 3.1. Erythrocyte membrane damage model The true stresses defined by the Abaqus built-in Mooney Rivlin model and the erythrocyte membrane damage model with the parameter values provided in Table 2 agreed well for the intact membrane, i.e. for β1 = 0, and for true strain between 0 and 1.1 (Figure 5 a). This agreement indicates that the developed VUMAT subroutine was accurately implemented. Table 2: Material parameter values to represent an intact erythrocyte membrane with the erythrocyte membrane damage model D1 (mm2/N) C10 (µN/mm2) C01 (µN/mm2) β1 12 152 15.2 0 When the damage parameter was changed from the intact state with β1 = 0 to the damaged state with β1 = 2.7, the in-plane true stress decreased for any given true strain (Figure 5 b). (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 28 For a true strain of 1.1, the in-plane true stress decreased from 0.0048 MPa for β1 = 0 to 0.0028 MPa for β1 = 0.49, 0.0016 MPa for β1 = 1, and 0.00025 MPa for β1 = 2.7. Figure 5: (a) True stress versus true strain in a shell element determined with the Abaqus built-in Mooney Rivlin model (‘Abaqus intact’) and the erythrocyte membrane damage model (‘Vumat intact’) for β1 = 0, and (b) true stress versus true strain determined with the erythrocyte membrane damage model showing the stress decrease with initiation and increase of erythrocyte membrane damage from β1 = 0 to 2.7. 3.2. Validation of the erythrocyte and invasion models The developed erythrocyte model was validated by comparing optical tweezer experimental data (Mills et al. 2004) with the numerical data obtained by simulating an optical tweezer experiment. When a force of 200 pN is applied diametrically, the dimension of the erythrocyte increases in the direction of the applied load, i.e. axial dimension (Figure 6 a) and decreases in the direction normal to the applied force, i.e. traverse dimension. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 29 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 30 Figure 6: a) The contour plots of axial and transverse displacement of the erythrocyte predicted with the FEM model at 0.0143 s, 0.209 s, 0.308 s, and 0.413 s of the optical tweezer simulation. b) Optical tweezer simulation data (axial and transverse erythrocyte finite element model diameters) fitted with optical tweezer experimental data for the axial and transverse diameter of a human erythrocyte from Mills et al. (2004). Contour plots of axial diameter (c) and maximum principal logarithmic strain (d) of the erythrocyte predicted with the finite element model of an optical tweezer test. The increase in displacement of the erythrocyte finite element model in the axial direction leads to an increase in the axial diameter of the erythrocyte. In contrast, the increase in the transverse displacements of the erythrocyte finite element model leads to the reduction of the transverse diameter of the erythrocyte (Figure 6 a). The numerical data shows that the transverse diameter of the erythrocyte model fit well with experimental data. However, the axial diameter of the erythrocyte model only fits well with experimental data when it is less than 0.0096 mm (Figure 6 b and c). For the erythrocyte model, the axial diameter of 0.0096 mm corresponds to the maximum principal erythrocyte membrane logarithmic strain of 1.81 (Figure 6 d). Beyond this point, the model fits experimental data in the axial direction with limited accuracy. The invasion model was validated by comparing the surface area of the erythrocyte with experimental erythrocyte areal deformation data obtained by tracking and segmenting the erythrocyte membrane during the invasion process (Geoghegan et al. 2021). The erythrocyte surface area numerically predicted without blood pressure fits well with the experimental data of the erythrocyte membrane from Geoghegan et al. (2021), also conducted without the application of blood pressure on the erythrocyte (Figure 7 a). The maximum error between the experimental areal data from Geoghegan et al. (2021) and the developed invasion model is 5.2%. Without blood pressure, the model predicts a decrease in the erythrocyte surface area consistent with experimental data from Geoghegan et al. (2021) study. In contrast, when a (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 31 blood pressure of 16 kPa is applied, the erythrocyte’s surface area increases during the late stages of the invasion process (Figure 7 a). The deformation of the erythrocyte model with a blood pressure of 16 kPa is entirely different from that without blood pressure (Figure 7 b and c). The largest maximum principal logarithmic strain in the erythrocyte membrane during the merozoite entry is 1.77 (Figure 7 d), which is less than the maximum principal logarithmic strain accuracy threshold of 1.81 established for the erythrocyte model. Figure 7: Comparison of numerically predicted (squares and diamonds) and experimental variation (circles) (Geoghegan et al. 2021) of erythrocyte surface area during merozoite invasion of an erythrocyte with membrane damage of β1 = 11 induced during the early invasion stage with τ = 0 to 0.1 s. Deformation of an erythrocyte at 100% indentation depth normalised over the length of the merozoite with (b) and without (c) blood pressure applied on the outer surface of the erythrocyte. d) Contour plot of the maximum principal logarithmic strain in the erythrocyte membrane with blood pressure at 100% normalised indentation depth. 3.3. Impact of erythrocyte morphology on merozoite invasiveness The impact of morphological variations of the erythrocyte on the invasiveness of the merozoite is determined by the maximum amount of energy required for successful invagination of the merozoite (case study 1). When the merozoite invades a convex region of the erythrocyte, the maximum strain energy predicted by the invasion model is 0.38 × 10-15 J, (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 32 and whereas when it invades the concave region of the erythrocyte, the maximum strain energy of 0.238 × 10-15 J is predicted (Figure 8 a). For invasion of a spherocyte, a maximum strain energy of 0.545 ×10-15 J is predicted, which is 43% and 129% larger than the energy required for the invagination of a merozoite in the convex and concave region of the erythrocyte, respectively. Figure 8: a) Required strain energy versus indentation depth predicted with the invasion model for invagination of a merozoite in the concave and convex region of the erythrocyte (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 33 and the spherocyte. b) Indentation force versus indentation depth for varying degrees and timing of erythrocyte membrane damage. 3.4. Impact of phosphorylation-induced damage in the erythrocyte membrane on the merozoite invasiveness The impact of malaria-induced erythrocyte membrane damage on the invasiveness of the merozoite (case study 2) is determined by the variation of the maximum indentation force required by the merozoite to invade the erythrocyte for increasing membrane damage represented by increasing β1 from 0 to 2.7 (Figure 8 b). For erythrocyte membrane damage, 𝑒𝛽1𝜏, limited to the early invasion stage (i.e. for τ = 0.1 s of the total simulation time of τ = 1.1 s), the maximum indentation forces are lower than for an equal amount of erythrocyte membrane damage induced for the entire simulation time. For example, the damage induced for β1 = 0.49 with τ = 0.1 s is the same as for β1 = 5.4 with τ = 1.1 s. Similarly, equal damage is obtained for β1 = 1 with τ = 1.1 s and β1 = 11 with τ = 0.1 s. The maximum indentation force is larger for an intact erythrocyte membrane (22.4 pN) than for an erythrocyte with membrane damage of β1 = 0.49 and τ = 1.1 s (17 pN). When the same damage amount is induced in the erythrocyte membrane during early invagination for τ = 0.1 s (with β1 = 5.4), the maximum indentation force is yet lower (15 pN). Similarly, the maximum indentation force of 11 pN for β1 = 11 with τ = 0.1 s is lower than the force of 12.5 pN for β1 = 1 with τ = 1.1 s. 3.5. Impact of local erythrocyte membrane damage on the global mechanical responses of the erythrocyte For compression of the entire erythrocyte (case study 3), the force-displacement curves do not differ between the normal erythrocyte and that with damage at location 1 and 2, (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 34 respectively, with β1 = 32 for τ = 0.1 s (Figure 9 a). This result shows that global compression is insensitive to local damage of the erythrocyte membrane due to phosphorylation. For local nanoindentation of an intact erythrocyte (case study 4), the maximum indentation force is 2.5-fold lower for damage induced at damage location 1 (i.e. the region of the indentation) (0.72 × 10-12 N) than at damage location 2 (1.78 ×10-12 N) (Figure 9 b). The latter is similar to the maximum indentation force for the intact erythrocyte membrane (1.68 ×10-12 N). An indentation displacement of 1.31 µm causes a maximum principal logarithmic strain in the erythrocyte membrane of 0.112 without blood pressure (Figure 9 c). For damage induced at damage locations 1 and 2, respectively, with β1 = 32 for 0.1 s, the maximum principal logarithmic strain is 0.589 (Figure 9 d) and 0.24 (Figure 9 e). The surface area of the intact and the damaged erythrocytes increases equally from 135 µm2 to 137 µm2, resulting in an areal strain of As,max = 0.014 (Figure 9 c to e). (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 35 Figure 9: Predicted force versus displacement for compression (a) and nanoindentation (b) of an intact erythrocyte and an erythrocyte damaged at damage locations 1 and 2, respectively. Maximum principal logarithmic strain for the intact erythrocyte membrane (c). Maximum principal logarithmic strain for the erythrocyte membrane with damage (β1 = 32) induced for 0.1 s at damage location 1 (d) and 2 (e) where A and As, max denotes the surface area and the areal ratio of the erythrocyte, respectively. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 36 4. Discussion Although there has been considerable progress in understanding which factors determine the merozoite’s invasiveness, most studies have not addressed the collective role played by the biophysical characteristics of the erythrocyte deformability in the invasion process. Cell shape, cytoplasmic viscosity, and membrane stability are the main determinants of erythrocyte deformability (Mohandas and Chasis 1993; Mohandas and Evans 1994). Experimental data on the mechanical properties of the human erythrocyte are widely available; however, it remains unknown how these physical properties influence the invasiveness of the merozoite. Hence, in the current study, an in silico approach was developed to investigate the role of erythrocyte morphology and merozoite-induced damage of the erythrocyte membrane in the invasiveness of the merozoite. 4.1. Erythrocyte membrane damage model To date, experimental investigations related to erythrocytes' invasion have predominantly focused on the role of parasite adhesins, signalling pathways, and the identity of binding receptors on the erythrocyte surface. Erythrocyte membrane damage mechanics associated with the invasion process has received limited attention (Koch et al. 2017). The erythrocyte membrane damage model developed in the current study utilises the hyperelastic Mooney Rivlin constitutive law, complemented with an exponential damage function to account for membrane remodelling due to phosphorylation. Despite limited knowledge of merozoite-induced damage, the model allowed the representation of various amounts of damage. One constraint of hyperelastic constitutive models is their instability when strain is inversely related to stress. One way to assess a model's stability is using Drucker’s criterion. However, this method has challenges since some hyperelastic models can be Drucker-unstable for small and large strains when subjected to different loading (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 37 conditions. Hence, common practice is to use a model with a known stability threshold or validity range. For example, the Mooney Rivlin model has a known validity range for equi- biaxial logarithmic strain up to 138% (Marckmann and Verron 2006). The stability threshold of the developed constitutive damage model was analytically determined using Drucker’s stability criteria. The decrease in the maximum stress represents stiffness reduction when damage is induced (Figure 5 b), indicating that the developed damage model accurately mimics the phosphorylation of the key erythrocyte membrane skeleton. The main limitation of the developed erythrocyte membrane damage model is the possible instability of the underlying Mooney Rivlin law at low strain. Furthermore, the second invariant of the left Cauchy tensor of the Mooney Rivlin law can make the developed model unstable in certain loading conditions. However, the decrease of the erythrocyte membrane stability threshold predicted by the developed model for increasing damage amount corresponds with reports that the phosphorylation of key protein elements of the erythrocyte skeleton leads to an unstable erythrocyte membrane. 4.2. Erythrocyte and invasion finite element models The developed erythrocyte model is a two-component model of the erythrocyte membrane and cytoplasm. The diametric dimensions of the erythrocyte predicted with the developed model for a simulated optical trap experiment were compared to experimental data from Mills et al. (2004) to validate the erythrocyte model. The predicted transverse erythrocyte diameters agree well with the experimental transverse diameters from Mills et al. (2004). However, the predicted axial diameter only fits the experimental data for an axial diameter smaller than 0.0096 mm, corresponding to the maximum principal logarithmic strain of the erythrocyte membrane of 1.81. Beyond this threshold value, the developed erythrocyte model has limited accuracy. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 38 The invasion model was validated by comparing the predicted surface area of the erythrocyte during the merozoite invasion with experimental data from Geoghegan et al. (2021). Since the erythrocyte was not exposed to blood pressure during the in vitro experiments, blood pressure was neglected in the invasion model for the validation simulation. The numerically predicted surface area excludes the region of the erythrocyte membrane that is in contact with the merozoite during wrapping. To allow direct comparison of the change in surface area during the invasion, the predicted data was normalised to the time of τ = 1.1, whereas the experimental data was normalised to the time of 49 s required for complete internalisation of the merozoite in the erythrocyte (Geoghegan et al. 2021). The predicted and experimental erythrocyte surface area agreed well (Figure 7 a). The validated model allows the prediction of the invasion process under in vivo conditions, e.g., blood pressure acting on the erythrocyte, which may be more challenging to achieve in in vitro experiments. The validation of the invasion model also allows the validation of the developed erythrocyte membrane damage model since the model predictions were in good agreement with experimental data for β1 = 11 and τ = 0.1 s. Identifying the damage parameters enables the identification of the correct invasion forces required by the merozoite to invade the erythrocyte. The invasion model indicated a maximum principal logarithmic strain in the erythrocyte membrane of 1.77. This value is below the accuracy threshold value of 1.81 of the erythrocyte model, indicating that the erythrocyte deformation during the invasion process is within the determined accuracy range of the erythrocyte model. 4.3. Impact of the erythrocyte morphology on the merozoite invasiveness The implication of morphological variations of the erythrocyte on the invasiveness of the merozoite was assessed by comparing the invasion energetics for (i) a discoid erythrocyte and a spherocyte and (ii) for the concave and convex membrane region of the discoid erythrocyte. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 39 The merozoite entry requires lower invasion energy for the discoid erythrocyte than the spherocyte. Reducing the surface area to volume ratio (S/V) increases the sphericity of the erythrocyte and leads to the formation of the spherocyte. The relatively low energy requirement indicates that the merozoite is more invasive when it invades a discoid-shaped erythrocyte than a spherocyte. An increase in sphericity corresponds to the increase in the energy required for the merozoite to invade the erythrocyte. The S/V of 1.44/m allows a healthy erythrocyte to undergo a large deformation of up to 230% of its original dimension. Reducing the healthy erythrocyte’s S/V by 14% forms a spherocyte with a surface area of 98.5 µm2 compared to the surface area of 135 µm2 of a healthy erythrocyte. The discoid erythrocyte shape provides an excess surface area of 36.5 µm2, i.e. 4.6-fold the surface area of 8.0 µm2 of a merozoite (Dasgupta et al. 2014), sufficient to facilitate the wrapping of the merozoite. The maximum strain energy predicted with the developed finite element invasion model corresponds to the total indentation work described by our analytical model (Msosa et al. 2023). The maximum strain energy of 38.0 × 10-17 J and 23.8 × 10-17 J predicted with the finite element invasion model for invasion in the convex and concave erythrocyte membrane region, respectively, is larger than the total indentation work of Ei = 1.40 × 10-17 J predicted by the analytical model for an areal strain of As,max = 51%. The higher strain energy predicted with the finite element invasion model may be due to the deformation of the erythrocyte cytoplasm, which is not considered in the analytical model. Erythrocytes with membrane protein abnormalities, such as hereditary spherocytosis, are generally spherical and less deformable than normal discoid erythrocytes. However, it is unknown whether these alterations may present the merozoite with a less ideal condition for invasion of erythrocytes by merozoites. Spherocytes have been found to have a low susceptibility to invasion by a merozoite. One of the reasons for low susceptibility is the genetic alteration of membrane proteins. However, from the current study, it has been (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 40 determined that shape alteration of the erythrocyte to spherical shape could be one of the contributing factors for the low susceptibility of spherocytes to infection by merozoites. 4.4. Impact of phosphorylation-induced damage in the erythrocyte membrane on merozoite invasiveness The impact of erythrocyte membrane damage on the invasiveness of the merozoite was studied by inducing local erythrocyte membrane damage. The amount of damage in the erythrocyte membrane model was regulated by varying damage parameter β1 between 0.49 and 2.7 such that β1 = 0.49 represented the minimum amount of damage and β1 = 2.7 represented the maximum amount of damage. The invasiveness of the merozoite was assessed by comparing the maximum indentation forces for each value of the damage parameter β1. The indentation force decreases with an increase in the amount of damage in the erythrocyte membrane model (Figure 9). This demonstrates that the invasiveness of the merozoite increases with the amount of damage. Merozoite-induced erythrocyte membrane damage has received limited attention, and erythrocyte membrane remodelling or damage stages are unknown. It is also unknown whether damage is induced only at the early invasion stage or throughout the invasion process. To validate the developed invasion model, erythrocyte membrane damage was induced at the beginning of the invasion process, i.e. at τ = 0.1 s with β1 = 11. The results suggest that erythrocyte membrane damage occurs during an early invasion process. The merozoite requires a greater force when damage is induced progressively throughout the invasion, i.e., for τ = 1.1 s, than at the beginning of the invasion process with τ = 0.1 s. These results demonstrate that the merozoite is more invasive when damage is induced during the early invasion stage (τ = 0.1 s) than progressively throughout (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 41 the invasion process. Hence, regulating the timing at which the merozoite induces erythrocyte membrane damage could be a potential target for antimalarial compounds. 4.5. Impact of local erythrocyte membrane on the global mechanical responses of the erythrocyte Compression simulations investigated the impact of local erythrocyte membrane damage on a global scale. The compression force does not differ for the intact and damaged erythrocyte, irrespective of the location of the phosphorylation damage (Figure 9 a). Hence, global compression of single erythrocytes cannot be successfully used to identify erythrocyte membrane damage. The simulations of nanoindentation of the erythrocyte in the central region indicate a discernible difference in the indentation force between an intact erythrocyte and an erythrocyte with membrane damage for damage in the central, concave region but not in the convex region of the discoid cell (Figure 9 e). This finding demonstrates that merozoite-induced local membrane damage may be detected with nanoindentation depending on the damage location, and further research is required. 5. Conclusions In this study, a finite element invasion model was developed and used to computationally quantify the mechanics of the invasion of a malaria merozoite into an erythrocyte and to investigate the impact of erythrocyte shape and membrane damage on the invasiveness of a malaria merozoite. The findings include the smallest force required for the malaria merozoite to invade a human erythrocyte successfully, i.e. 11 pN. The invasiveness of the merozoite decreases with an increase in the sphericity of the erythrocyte, which is associated with genetic disorders such as hereditary spherocytosis. An increase in phosphorylation-induced membrane damage in the erythrocyte increases the invasiveness of the malaria merozoite, as may be expected. It was further found that the malaria merozoite is more invasive when (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 42 erythrocyte membrane damage induced by phosphorylation is limited to an early invasion stage compared to the entire invasion stage. The findings on the invasion mechanics can guide future experimental studies to assess the invasiveness of the merozoite. The results from the nanoindentation simulations indicate the suitability of nanoindentation as an additional experimental technique to determine erythrocyte membrane damage in the context of invasion-blocking anti-malaria drugs. The developed computational models of the human erythrocyte and merozoite invasion can be adapted to study other parasite invasion processes. Funding This research was supported financially by the National Research Foundation of South Africa (grants CPRR14071676206 and IFR14011761118 to TF) and the South African Medical Research Council (grant SIR328148 to TF), and grants from the World Bank to the University of Malawi. The funders had no role in study design, data collection and analysis, the decision to publish, or the preparation of the manuscript. Any opinions, findings, conclusions, or recommendations expressed in this publication are those of the authors and do not necessarily represent the official views of the funding agencies. Conflicts of Interest The authors declare no conflict of interest. Data availability Software used and data supporting the results presented in this article are available on the University of Cape Town’s institutional data repository (ZivaHub) under https://doi.org/10.25375/uct.28263767 as Msosa C, Abdalrahman T, Franz T. Software code and data for “In silico analysis of the invasion mechanics and invasiveness of the plasmodium falciparum merozoite", Cape Town, ZivaHub, 2025, DOI 10.25375/uct.28263767. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 43 Reviewers may access the data record with this private link https://figshare.com/s/a5c508ec78a4af853be6. This statement will be omitted in the accepted manuscript. CRediT author contributions CM: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – Original Draft, and Writing - Review & Editing TA: Conceptualization, Methodology, Project administration, Supervision, and Writing - Review & Editing TF: Conceptualisation, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Visualization, and Writing - Review & Editing

Abdalrahman T, Franz T (2017) Analytical modeling of the mechanics of early invasion of a merozoite into a human erythrocyte Journal of biological physics 43:471-479 doi:10.1007/s10867-017-9463-6 Dasgupta S et al. (2014) Membrane-Wrapping Contributions to Malaria Parasite Invasion of the Human Erythrocyte Biophys J 107:43-54 doi:10.1016/j.bpj.2014.05.024 Eber S, Lux SE (2004) Hereditary spherocytosis--defects in proteins that connect the membrane skeleton to the lipid bilayer Semin Hematol 41:118-141 doi:10.1053/j.seminhematol.2004.01.002 Flannery EL, Chatterjee AK, Winzeler EA (2013) Antimalarial Drug Discovery: Approaches and Progress towards New Medicines Nat Rev Microbiol 11:849-862 doi:10.1038/nrmicro3138 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 44 Geoghegan ND et al. (2021) 4D analysis of malaria parasite invasion offers insights into erythrocyte membrane remodeling and parasitophorous vacuole formation Nature Communications 12:3620 doi:10.1038/s41467-021-23626-7 Gilson PR, Crabb BS (2009) Morphology and kinetics of the three distinct phases of red blood cell invasion by Plasmodium falciparum merozoites Int J Parasitol 39:91-96 doi:10.1016/j.ijpara.2008.09.007 Gingold RA, Monaghan JJ (1977) Smoothed particle hydrodynamics: theory and application to non-spherical stars Mon Not R Astron Soc 181:375-389 doi:10.1093/mnras/181.3.375 Hochmuth RM, Mohandas N, Blackshear PL (1973) Measurement of the Elastic Modulus for Red Cell Membrane Using a Fluid Mechanical Technique Biophys J 13:747-762 Koch M et al. (2017) Plasmodium falciparum erythrocyte-binding antigen 175 triggers a biophysical change in the red blood cell that facilitates invasion Proceedings of the National Academy of Sciences of the United States of America 114:4225-4230 doi:10.1073/pnas.1620843114 Li W (2016) Damage Models for Soft Tissues: A Survey Journal of Medical and Biological Engineering 36:285-307 doi:10.1007/s40846-016-0132-1 Li X, Peng Z, Lei H, Dao M, Karniadakis GE (2014) Probing red blood cell mechanics, rheology and dynamics with a two-component multi-scale model Philos Trans A Math Phys Eng Sci 372 doi:10.1098/rsta.2013.0389 Li Y, Lu L, Li J (2016) Topological Structures and Membrane Nanostructures of Erythrocytes after Splenectomy in Hereditary Spherocytosis Patients via Atomic Force Microscopy Cell Biochem Biophys 74:365-371 doi:10.1007/s12013-016-0755- 4 Marckmann G, Verron E (2006) Comparison of Hyperelastic Models for Rubber-Like

Rubber Chem Technol 79:835-858 doi:10.5254/1.3547969 %J Rubber Chemistry and Technology (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 45 Mills JP, Qie L, Dao M, Lim CT, Suresh S (2004) Nonlinear elastic and viscoelastic deformation of the human red blood cell with optical tweezers Mechanics and Chemistry of Biosystems 1:169-180 Mohandas N, Chasis JA (1993) Red blood cell deformability, membrane material properties and shape: regulation by transmembrane, skeletal and cytosolic proteins and lipids Semin Hematol 30:171-192 Mohandas N, Evans E (1994) Mechanical properties of the red cell membrane in relation to molecular structure and genetic defects Annu Rev Biophys Biomol Struct 23:787-818 doi:10.1146/annurev.bb.23.060194.004035 Monaghan JJ (1988) An introduction to SPH Comput Phys Commun 48:89-96 doi:https://doi.org/10.1016/0010-4655(88)90026-4 Msosa C, Abdalrahman T, Franz T (2023) An analytical model describing the mechanics of erythrocyte membrane wrapping during active invasion of a plasmodium falciparum merozoite J Mech Behav Biomed Mater 140:105685 doi:https://doi.org/10.1016/j.jmbbm.2023.105685 Ogden RW (1978) Nearly isochoric elastic deformations: Application to rubberlike solids J Mech Phys Solids 26:37-57 doi:https://doi.org/10.1016/0022-5096(78)90012-1 Pinder J, Fowler R, Bannister L, Dluzewski A, Mitchell G (2000) Motile Systems in Malaria Merozoites: How is the Red Blood Cell Invaded? Parasitol Today 16:240-245 doi:10.1016/S0169-4758(00)01664-1 Preiser P, Kaviratne M, Khan S, Bannister L, Jarra W (2000) The apical organelles of malaria merozoites: host cell selection, invasion, host immunity and immune evasion Microbes and Infection 2:1461-1477 doi:https://doi.org/10.1016/S1286- 4579(00)01301-0 Riglar DT et al. (2011) Super-resolution dissection of coordinated events during malaria parasite invasion of the human erythrocyte Cell Host Microbe 9:9-20 doi:10.1016/j.chom.2010.12.003 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint 46 Song H, Liu Y, Zhang B, Tian K, Zhu P, Lu H, Tang Q (2017) Study of in vitro RBCs membrane elasticity with AOD scanning optical tweezers Biomedical Optics Express 8:384-394 doi:10.1364/BOE.8.000384 Systemes D (2015) Abaqus 6.14 documentation–theory guide Providence, RI Vieira AC, Marques AT, Guedes RM, Tita V (2011) Material model proposal for biodegradable materials Procedia Engineering 10:1597-1602 doi:10.1016/j.proeng.2011.04.267 Violeau D, Rogers BD (2016) Smoothed particle hydrodynamics (SPH) for free-surface flows: past, present and future Journal of Hydraulic Research 54:1-26 doi:10.1080/00221686.2015.1119209 Wang C, Wang Y, Peng C, Meng X (2016) Smoothed particle hydrodynamics simulation of water-soil mixture flows Journal of Hydraulic Engineering 142:04016032 Weiss JA (1994) A Constitutive Model and Finite Element Representation for Transversely Isotropic Soft Tissues. University of Utah Ye T, Phan-Thien N, Lim CT (2016) Particle-based simulations of red blood cells—A review Journal of Biomechanics 49:2255-2266 Zhang Z, Zhang X (2011) Mechanical behavior of the erythrocyte in microvessel stenosis Science China Life Sciences 54:450-458 doi:10.1007/s11427-011-4152-3 Zuccala ES, Satchwell TJ, Angrisano F, Tan YH, Wilson MC, Heesom KJ, Baum J (2016) Quantitative phospho-proteomics reveals the Plasmodium merozoite triggers pre- invasion host kinase modification of the red cell cytoskeleton Sci Rep 6:19766 doi:10.1038/srep19766 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted June 27, 2025. ; https://doi.org/10.1101/2025.06.26.661885doi: bioRxiv preprint