In-Flow dynamics of an area-difference-energy spring-particle red blood cell model on non-uniform grids

Abstract In this paper the area-difference-energy spring-particle (ADE-SP) red blood cell (RBC) structural model developed by Chen and Boyle is coupled with a lattice Boltzmann flux solver to simulate RBC dynamics. The novel ADE-SP model accounts for bending resistance due to the membrane area difference of RBCs while the lattice Boltzmann flux solver offers reduced computational runtimes through GPU parallelisation and enabling the employment of non-uniform meshes. This coupled model is used to simulate RBC dynamics and predictions are compared with existing experimental measurements. The simulations successfully predict tumbling, tank-treading, swinging and intermittent behaviour of an RBC in shear flow, and demonstrate the capability of the model in capturing in-flow RBC behaviours.


Introduction
In recent years, there have been massive developments in the computational modelling of red blood cell (RBC) dynamics. A variety of RBC dynamic behaviours have been predicted including tumbling(TB), tank-treading (TT), and swinging dynamics in simple shear flows (Noguchi and Gompper 2005a), slipper and parachute shapes encountered in capillary flow (Fedosov et al. 2014), and multi-cell behaviours such as the formation of rouleaux aggregates (Fedosov et al. 2014) and cell free layers (Fedosov et al. 2014). The two main approaches to the modelling of the structural behaviour of a RBC are the continuum approach and the discrete particle approach. What they have in common is the treatment of the membrane as a surface. The composite membrane of a RBC has a thickness of approximately 9-10 nm; this is orders of magnitude lower than the length scale of microns for a RBC. This allows the membrane to be treated as a two-dimensional (2D) surface embedded in three-dimensional (3D) space. The continuum approach predicts the structural mechanics of a solid by applying constitutive laws which define energy-strain relationships. Many laws have been suggested to model the in-plane mechanics of biological membranes including the Skalak membrane (Skalak et al. 1973) and the Mooney-Rivlin membrane (Li et al. 1988). Spring-network models can also be used to model the solid mechanics of a surface including RBC membranes. This approach involves discretising the surface into particles which are interconnected with springs. The mechanical properties of these springs then directly influence the behaviour of the model. The most common type of spring used to model biological membranes are worm like chains (WLC) to model the in-plane shear resistance of the membrane (Discher et al. 1998). As with the continuum approach, there are also differing approaches to modelling the out-of-plane bending resistance of the RBC. These include the spontaneous curvature model (Discher et al. 1998) and the area-difference-energy (ADE) approach (Chen and Boyle 2017), with the latter accounting for non-local bending resistance due to the area difference between the outer and inner layers of the lipid bilayer.
The earliest attempts at modelling RBC membrane mechanics looked at spring-networks at the spectrin level (Discher et al. 1998). This resulted in accurate predictions that compared favourably with experimental data and predictions of continuum structural solvers.
However this approach is extremely computationally expensive with upwards of 20,000 vertices required. More recently coarse grained spring-particle (CG-SP) models have been developed where a systematic coarse-graining procedure enables a reduction in the number of degrees of freedom in the RBC model by two orders of magnitude (Noguchi and Gompper 2005b;Pivkin and Karniadakis). Similar to the continuum models, the CG-SP approach has successfully modelled the RBC mechanics in two-phase flow. An advantage of the CG-SP methods over continuum methods is the ability to account for membrane viscosity, thermal fluctuations and diseased RBC dynamics (Fedosov et al. 2011).
In this work the CG-SP model with the ADE bending resistance of Chen and Boyle (2017) is adopted and coupled with a Navier-Stokes solver to simulate RBC dynamics in fluid flow. The ADE model is necessary to account for the bending resistance due to the membrane area difference (MAD) of the RBC in a CG-SP model. A goal of this work is to show that RBC dynamics can be captured while accounting for the MAD of the RBC. In addition, a significant challenge with modelling RBC dynamics is the limitation on the physical time that can be modelled per timestep in a numerical model. Ways to reduce this constraint include using more computational processing power and reducing the number of fluid cells via the employment of non-uniform/unstructured grids. To fully utilise modern computational processing power, a highly parallelisable numerical method is required. Incorporating the lattice Boltzmann flux solver into RBC dynamical modelling brings a highly parallelisable numerical method that can scale on GPUs while successfully employing non-uniform grids. In the rest of this paper the governing equations of the ADE-CG-SP solver are introduced, the simulations performed with the model are described, and comparisons are made of predictions with experimental measurements from the literature.

Governing equations
In this section a brief overview of the governing equations used in this work is given. A detailed description of the discretised version of the governing equations is provided in the Supplementary Material. The suspending fluid was modelled using the incompressible isothermal Navier-Stokes equations and predictions of the fluid flow were made using a lattice Boltzmann flux solver (LBFS) (Wang et al. 2015;Walsh and Boyle 2020). The LBFS uses a local solution of the lattice Boltzmann equation to calculate the fluxes between cell volumes. The benefit of the LBFS is that it doesn't involve the expensive solution of the Poisson equation for pressure correction while being applicable to nonuniform and unstructured meshes. Due to their computational efficiency, non-uniform Cartesian grids were used to discretise the computational domain.
An immersed boundary method (IBM) was used to model the fluid-structure interaction between the fluid and the RBC. Jang and Lee (2017) used a reproducing polynomial particle method (RPPM) to calculate the velocity interpolation and force spreading calculations of the IBM on non-uniform grids. This method generates a linear system of equations to ensure the conservation of moments in the system. However it doesn't involve any rescaling or arbitrary parameters. One massive advantage of the RPPM is that the moment matrix is in the form of a Vandermonde matrix and allows the linear system to be solved explicitly. For this reason the method of Jang and Lee (2017) was adopted in this work.
The RBC consists of a surface membrane which surrounds an interior fluid. This fluid, known as cytosol, is an incompressible fluid which ensures that volume is conserved. The surface membrane consists of two main parts which confer structural strength to the RBC: the plasma membrane (PM) and the cytoskeleton. The PM provides structural resistance to surface compressibility and out-of-plane bending, while the cytoskeleton provides resistance to in-plane shearing. Each of these mechanical properties provides a corresponding contribution to the total RBC Helmholtz energy, E RBC , defined as: where E V and E A are the energies due to the volume and area constraints respectively, E B is the bending energy, and E S is the spring potential energy (Discher et al. 1998). Key to the ADE-CG-SP RBC model is that the bending energy term has two contributions. The first is the Helfrich term which is due to local bending resistance. The second is the ADE term which is due to the contribution from the membrane area difference (MAD). Including this ADE term allowed Chen and Boyle (2017) to accurately model the complex stomatocyte-dicsocyte-echinocyte transformation. In the ADE-CG-SP RBC model, springs are modified WLC, with the spring potential energy term defined as in Chen and Boyle (2017) as: where k B is the Boltzmann constant, T K is the absolute temperature in Kelvin, L P is the persistent length of the WLC springs, s is the spring index, N s is the total number of springs in the RBC mesh, c 1 and c 2 are constants, and k s is the ratio of the instantaneous length,L s , of the spring to the contour length, L c, s , of spring s, i.e., The persistent length is the length at which the WLC transitions from stiff to flexible behaviour and the contour length is the maximum permissable length of a WLC. The conservative force which governs the structural behaviour of the RBC is defined as: where R is the vector of particle coordinates in the RBC. Note that the force is in the direction that will cause the greatest reduction in free energy. Equation (4) also shows that the conservative force is dependent on the energy contribution from each of the four mechanical properties, i.e., area, volume, bending and shear moduli. The values of the mechancial properties used are given in Table 1. To account for the influence of membrane viscosity, the approach of Ehi-Egharevba (date unknown) is applied and is described in the Supplementary Material.

Methods
In this work the blood flow model incorporating a LBFS, a non-uniform-grid IBM and an ADE-SP RBC model was validated against experimental measurements for individual RBCs. These experiments included:

Optical Tweezers Test In-Plane Shear Flow -Wheel Configuration
Out-of-Plane Shear Flow -Tumbling, Tank-Treading and Swinging In this work, all three experiments were first performed using the stress free RBC as used by Chen and Boyle (2017). It was only during the tank-treading simulation that it was discovered that a prestressed RBC was required to recover physiologically realistic results. The optical tweezers and wheel configuration were then repeated with the pre-stressed configuration to confirm its suitability for modelling RBC mechanics. In this paper the RBC was prestressed by prescribing a persistent length to the WLC spring as done by Pivkin and Karniadakis (2008).

Common Mesh types
Each of the three benchmark flow problems had a common approach employed to meshing the suspending fluid. A RBC was immersed in a cubic fluid domain with varying edge length. The fluid was meshed using a non-uniform Cartesian mesh with a varying number of cells. A plateau function was used to generate the distribution of the cell size in the x, y and z directions. Please refer to the Supplementary Material for details of the plateau functions.

RMS error
The root mean square (RMS) of the error between experimental measurements and numerical predictions is defined as: and is used as a metric to quantify the accuracy of the numerical predictions.

Introduction
The optical tweezers experiment was performed by Mills et al. (2004) to investigate the mechanical properties of RBCs in response to stretching. The RBCs in this experiment were procured from healthy adults using a lancet device and isolated using a centrifugation process. Two silica micro beads 4.12 mm in diameter were then binded to diametrically opposite points on the RBC. This took place in a suspending medium of phosphatebuffered-saline (PBS) with PH 7.4 In this experiment the optical tweezers comprised a single-bead gradient optical trap. This involved the use of a laser to optically trap one of the silica beads. The second silica bead was then adhered to the glass slide on the stage. To apply a force the stage was then translated and a resisting force up to a maximum of 193 ± 20 pN was provided by the laser (see Figure 1 for an illustration of the experiment).

Problem set-up
The following approach was taken to numerically simulate this experiment. The computational domain used was cubic with an edge length of 24 mm. The mesh was generated as described in Section 3.1 with 40 cells used along each of the Cartesian axes with plateau function values of a ¼ 24 and b ¼ 0.
2. An illustration of this mesh is provided in Figure 2. This results in cell heights of 0.3 mm in the vicinity of the RBC-fluid interface while using 64,000 cells in the mesh. An equivalent uniform mesh would require 512,000 cells. The RBC mesh had 2,562 nodes and its initial shape was taken as a relaxed biconcave discocyte described in the Supplementary Material and was positioned at the centre of the fluid domain as per Figure 2. Chen and Boyle (2017) showed that there was very little variation in results for the optical tweezers test with mesh density once a minimum of 289 nodes was used. PBS is considered to have similar mechanical properties to water and a viscosity ratio k l ¼ l cytosol =l externalfluid %5 was used in this test case. The fluid had initial conditions of V ¼ 0 and q fluid ¼ 1000 kg m À3 : Neumann boundary conditions specifying zero gradient across the boundary were applied to all faces for both velocity and density. In Mills et al. (2004) experiment a range of stretching forces, F Stretching , from 0 to 193 pN were applied to the RBC via the attached silica beads. The contact diameter D c between the silica beads and the RBC was 2 mm (Dao et al. 2003). F Stretching was then distributed over 4% of the nodes on the RBC mesh to reflect the contact diameter of 2 mm. This is 2% of the nodes with the largest values of x-coordinate and 2% of the nodes with the smallest values of x-coordinate. The absolute value of the force at each node F Stretching, i ¼ F Stretching 0:02N i but orientated in opposite directions for maximum and minimum x-coordinate nodes (see Figure 3).
After F Stretching was applied, the simulation was run until equilibrium was reached. Then the axial diam- The simulation was performed twice, first using the RBC parameters for a stress free biconcave RBC from Table 1. This was done to show that the present work could recreate the results of Chen and Boyle (2017) when the RBC was immersed in a suspending fluid. The simulation was then performed a second time with the parameters for a pre-stressed biconcave RBC from Table 1.

Results
The optical tweezers experiment was simulated for F Stretching ¼ 0À200 pN with simulations performed at 10 pN intervals for a total of 20 simulations. The D A and D T at static equilibrium are plotted against F Stretching in Figure 4; as can be seen there is excellent agreement between the numerical predictions and experimental measurements for the stress free biconcave RBC. In Figure 5, the final experimental and   steady-state numerical predictions are plotted for a variety of forces. Again there is excellent agreement between both sets of results. It is apparent that as the force increases that D A increases and D T decreases. The results for a pre-stressed biconcave RBC are also plotted in Figure 4. As shown by the RMS error for stress-free and pre-stressed predictions in Table 2, these results are not as accurate as the stress free biconcave RBC for D A but are well within the error bars of the measurements of Mills et al. (2004). It is noticeable that the pre-stressed biconcave RBC D T predictions have better agreement with measurements.
5. Deformation in "wheel" configuration 5.1. Introduction Yao et al. (2001) investigated the efficacy of Low Viscosity Ektacytometry (LVE) by measuring the deformation of RBCs in a "wheel" configuration in a shear flow of low viscosity. The "wheel" configuration is when both major axes of the biconcave RBC are in the same plane as the shear flow applied (see Figure  6). The measurements from LVE were compared to those produced by a flow chamber and imaging system by using a deformation index, DI, which is defined as: where D max and D min are the lengths of the major and minor axes of the elliptical RBC. The imaging system involved taking 20 exposures of digital photographs of a RBC suspension in shear flow and measuring the DI directly from the photo of the RBCs in the "wheel" configuration. In comparison LVE involves measuring the DI by passing a laser through a RBC suspension in shear flow and measuring the elliptical diffractive pattern on a phototube. Experiments were performed for shear rates in the range 0s À1 <G<120s À1 : The suspending medium used was PBS with a viscosity ratio k l %8:48: For both approaches a linear relationship between DI and G was discovered: In this validation problem these linear relationships are used to measure the accuracy of the ADE-SP RBC model.

Problem Set-Up
A RBC was immersed in a cubic fluid domain with an edge length of 24 mm. The fluid domain was meshed using a non-uniform structured mesh with 8,000 cells. In this flow problem 20 cells were used along each of the Cartesian axes with plateau function values of a ¼ 8 and b ¼ 0.2. An illustration of this mesh is provided in Figure 6. This results in cell heights of 1.0 mm in the vicinity of the RBC-fluid interface while using 8,000 cells in the mesh. The initial RBC mesh had 2,562 nodes. Its initial shape was a biconcave RBC and was positioned at the centre of the fluid domain as per Figure 6. PBS is considered to have similar mechanical properties to water and a viscosity ratio k l %8:48 was used in this test case. The fluid had initial conditions of V ¼ 0 and q fluid ¼ 1000 kg m À3 : Dirichlet boundary conditions specifying a velocity equivalent to the shear rate G were applied to all faces. Neumann boundary conditions specifying zero gradient across the boundary were applied to all faces for density. In this work simulations were run in the range 0s À1 <G<120s À1 , with simulations performed at 20 s À1 intervals. The RBC had the material properties described in Table 1; however a sensitivity analysis was performed on the shear  Simulations were performed in the range of 2lN=m<K S <6:1lN=m: Simulations were also performed with the parameters for a pre-stressed biconcave RBC from Table 1.

Results
An initial simulation was run for three periods of revolution and it was found that 1.25 periods of revolution were sufficient to find a constant DI. The remaining simulations were run for 1.25 periods of revolution of the RBC. For each simulation the DI was calculated by extracting the x and y coordinates of the RBC mesh and finding the minimum sized rectangle that bounded these coordinates. The resulting DI is compared to the linear relationships of Yao et al. (2001) in Equations (7) to (8) in Figure 7. As can be seen the predictions show an approximate linear relationship between the DI and G for all values of K S . This is in agreement with the experimental work done by Yao et al. (2001). The RMS error of the numerical predictions defined in Equation (5) is shown in Table 3. The results with K S ¼ 4 lN=m with an initial stress free shape show a DI which is slightly lower than that of the experimental work. As seen from the sensitivity analysis, a value of K S ¼ 3 lN=m for a stress free shape has the lowest RMS error with a value of 0.137. When a pre-stressed initial shape was used with a value of K S ¼ 6:1lN=m, the RMS error is equal to 0.438. This is much larger  than K S ¼ 3 but is more accurate than the stress free results of K S ¼ 6:1lN=m which has an RMS error equal to 0.764.

Introduction
The dynamic behaviour of RBCs in shear flow, where the two major axes of the RBC are out of plane with the shear flow applied, has been the subject of many experimental studies. The observations from these studies have shown three primary dynamical modes: tumbling (TB) at low shear rates, tank treading (TT) at high shear rates, and an intermediate regime when transitioning from TB to TT (Fischer et al. 1978;Tran-Son-Tay et al. 1984). TB is where the pitch angle of the RBC, see Figure 8, undergoes continuous full 360 degree revolutions. TT is where the pitch angle remains approximately constant but the surface of the RBC undergoes a continuous revolution around the centroid of the RBC. The transition from TB to TT is due to the existence of" shape memory" in the RBC; this is the existence of a minimum energy when the RBC is in static equilibrium (Fischer 2004). The transition to TT occurs when the shear stress applied by a fluid exceeds the maximum elastic energy storage of the RBC (Skotheim and Secomb 2007). This transition point is a function of the shear rate G, viscosity ratio k l , and the elastic properties of the RBC. In experiments to date, isolated RBCs have been subjected to shear flow whilst suspended in various solutions of dextran. The solutions of dextran have varying viscosities and images were recorded of the dynamic behaviour of the RBCs. Numerical simulations of a RBC in shear flow were performed in this work and compared to the experimental work of Abkarian et al. (2007) and the numerical work of Pivkin and Karniadakis (2008) and Reasor et al. (2012).

Problem Set-Up
A RBC was immersed in a cubic fluid domain. The domain size was reduced compared to the previous flow problems due to the increased computational cost of this problem. Various domain sizes were trialled and an edge length of 12 mm was found to be the smallest edge length that did not adversely impact on the simulations performed. The fluid domain was meshed using a Cartesian non-uniform mesh with 27,000 cells. In this flow problem 30 cells were used along each of the Cartesian axes with plateau function values of a ¼ 10 and b ¼ 0.
2. An illustration of this mesh is provided in Figure 9. The initial RBC mesh had 2562 nodes and its initial shape was taken as the pre-stressed biconcave RBC described in the Supplementary Material and was positioned at the centre of the domain as per Figure 9. Dextran was considered to have similar mechanical properties to water and the fluid had initial conditions of V ¼ 0 and V ¼ 0: In this work, the approach of the experimental work of Abkarian et al. (2007) and the numerical work of Pivkin and Karniadakis (2008) and Reasor et al. (2012) was adopted where a suspending fluid with a viscosity of 22 mPa s was used. This gave a viscosity ratio of k l ¼ 0:27: Dirichlet boundary conditions specifying a velocity equivalent to the shear rate G were applied to all faces for all velocity components. Neumann boundary conditions specifying zero gradient across the boundary were applied to all faces for density. In this work, simulations were run in the range 0s À1 <G<7:5s À1 : The RBC had the material properties described in Table 1.

Results
The first benchmark parameter identified was the period of rotation of the RBC under varying shear rates. It was also identified if the RBC was TB or TT Table 3. Root mean square of the error of numerical predictions of the wheel experiment versus experimental measurements of Yao et al. (2001). at each shear rate. The results are shown in Figure 10. As can be seen there is very good agreement between the present predictions and the predictions of Pivkin and Karniadakis (2008) and Reasor et al. (2012). The present work correctly predicts that lowering the shear rate results in longer periods of oscillation of the RBC. The present work also predicts a transition between TB and TT in the range of 1s À1 <G<1:4s À1 : Reasor et al. (2012) predicted a range of 1s À1 <G<1:5s À1 with Pivkin and Karniadakis (2008) predicting a range of 1:15s À1 <G<1:35s À1 : Snapshots of a RBC TB (G ¼ 0:5s À1 ) and a RBC TT (G ¼ 1:6s À1 ) are shown in Figures 11 and 12 respectively. As is shown, the TB RBC fully rotates around the z Cartesian axis. In comparison, the TT RBC slightly oscillates around a fixed pitch angle from the xz Cartesian plane. This behaviour is referred to as swinging and the amplitude of this angle has been experimentally measured by Abkarian et al. (2007) for the case of G ¼ 1:8s À1 and compared with predictions in Figure 13. It can be seen that there is very good agreement between the present work and the experimental measurements of Abkarian et al. (2007) for the amplitude of the pitch angle. The initial period of the results should be discounted as the RBC starts from steady state and initial perturbations are removed. Abkarian et al. (2007) also discovered a logarithimic relationship between the amplitude of the swinging angle and the shear rate. The results from the present work are shown in Figure 14. The linear curve fitted to the results has a slope of À0.983 which is in excellent agreement with Abkarian's slope of À1. Finally, intermittent behaviour was also observed in the present work for G ¼ 1:2s À1 and is shown in Figure 15. As can be seen, the RBC alternates between TB and TT behaviour. Abkarian et al. (2007) Pivkin and Karniadakis (2008) and Reasor et al. (2012). The mode of the rotation of the RBC is also denoted. behaviour at a much higher value of G ¼ 1:526s À1 : This was not observed in the current work. The numerical works of Reasor et al. (2012) and Pivkin and Karniadakis (2008) did not show results for intermittent behaviour of a RBC but claim the transition/ intermittency zone is between 1s À1 <G<1:5s À1 and 1:15s À1 <G<1:35s À1 respectively.

Flow problem runtimes
The runtimes for the flow problems are shown in Table 4. These runtimes confirm that the allowable time step is the primary factor in determining the computational effort required. The timestep Dt is of the order of ms whereas the flow problem requires time scales in the order of seconds to be simulated. This suggests that an alternate approach to time integration that allows larger time steps would offer large computational savings. It is also noticeable that the time per iteration is higher for the wheel problem than the TT problem even though the wheel problem has a lower N cells . The increase in computational effort per iteration in the wheel problem is due to the VI-FS steps during the IBM process. This demonstrates the complex relationship between computational complexity and runtimes when GPUs are used. Depending on the orientation and position of the RBC this could lead to uncoallesced memory accesses and code divergence.

Summary
This paper has shown that the ADE-SP LBFS blood flow model can capture RBC dynamics while accounting for the additional bending resistance due to the MAD of the RBC. This will enable further research into the behaviour of more complicated RBC morphologies in shear flow. It should enable the experimentally measured TB-TT transitions identified by Fischer and Korzeniewski (2013) for echinocytes and stomatocytes to be numerically predicted. The results presented in this paper provide great confidence that this can be achieved. The optical tweezers numerical experiment performed produced equivalent results to Chen and Boyle (2017) with the exact same RBC parameters. However when the ADE-SP LBFS blood flow model was validated against further experimental results, it was found that a pre-stressed biconcave Figure 11. Snapshots of a RBC tumbling over a period where G ¼ 0:5s À1 and k l ¼ 0:27: The normalised velocity magnitude of the suspending fluid is shown. The node index of the mesh is also shown to demonstrate the rotation of the RBC.
geometry was required to simulate RBC dynamics with physiologically realistic values. A pre-stressed geometry using the values in Table 1 produces very accurate results for RBC dynamics in shear flow that agree with experimental evidence in the literature. However this produced more inaccurate results than the stress free configuration when the pre-stressed    Abkarian et al. (2007) between the variables is also plotted.
configuration was retrospectively used to simulate the optical tweezers and wheel deformation experiment. It resulted in larger deformations in the optical tweezers experiment and smaller deformations in the wheel experiment. This may be down to the spring model used. The modified WLC model proposed by Chen is stiffer than alternative spring models used in the literature. Fedosov (2010) used a WLC-POW model whereas Pivkin and Karniadakis (2008) and Reasor et al. (2012) used a WLC spring mixed with a hydrostatic elastic energy term. Another consideration is the choice of model for membrane viscosity. Using the viscosity model of Ehi-Egharevba (date unknown) does not impact on RBC dynamics in blood flow scenarios. The work of Fedosov et al. (2014) showed that membrane viscosity has significant impact on the RBC tank treading frequency and also swinging angle amplitudes. Membrane viscosity is captured in this work by assigning the membrane viscosity to the SPH particles in the fluid colocated with the membrane. This would suggest that to capture membrane viscosity in the current work, that the fluid volume's viscosity must be updated to reflect the viscosity of the membrane. On the other hand Pivkin and Karniadakis (2008) and Reasor et al. (2012) do not supply results for the amplitude of swinging angles in their work. It may be that if these models were adopted that they do not accurately capture the pitch orientations or swinging angles that are observed in experimental studies. The current work also shows that the LBFS is a suitable solver for use in blood flow simulations. This solver could potentially offer a very cost effective approach to modelling interaction of RBCs with complex geometries of biomedical devices due to its scalable parallelisation and mesh refinement capabilities. Finally, the runtimes of the simulations performed confirm that the allowable time step is the primary factor in determining the computational effort required. The timestep Dt is of the order of ms whereas the flow problem requires time scales in the order of seconds to be simulated. This would suggest that an alternate approach to time integration that allows larger time steps would offer large computational savings. Future work required to improve this model include the inclusion of cell-tocell interaction forces. This would allow the inclusion of other cells in the computational domain e.g., white blood cells and platelets. These cells could be modelled using the existing spring-network model or with bespoke models. Additional future work required to improve this model includes the ability to handle multi-zone fluid domains. This would enable an optimal approach to solving bifurcation, bended artery and constriction problems as found in Lakshmanan et al. (2019) and Javadzadegan et al. (2016). The ability to handle unstructured meshes currently exists in the model but was not used in this work to optimise runtimes.