Macrophage Anti-inflammatory Behaviour in a Multiphase Model of Atherosclerotic Plaque Development

Atherosclerosis is an inflammatory disease characterised by the formation of plaques, which are deposits of lipids and cholesterol-laden macrophages that form in the artery wall. The inflammation is often non-resolving, due in large part to changes in normal macrophage anti-inflammatory behaviour that are induced by the toxic plaque microenvironment. These changes include higher death rates, defective efferocytic uptake of dead cells, and reduced rates of emigration. We develop a free boundary multiphase model for early atherosclerotic plaques, and we use it to investigate the effects of impaired macrophage anti-inflammatory behaviour on plaque structure and growth. We find that high rates of cell death relative to efferocytic uptake results in a plaque populated mostly by dead cells. We also find that emigration can potentially slow or halt plaque growth by allowing material to exit the plaque, but this is contingent on the availability of live macrophage foam cells in the deep plaque. Finally, we introduce an additional bead species to model macrophage tagging via microspheres, and we use the extended model to explore how high rates of cell death and low rates of efferocytosis and emigration prevent the clearance of macrophages from the plaque.


Introduction
Atherosclerosis is a chronic inflammatory disorder characterised by the retention of lipoproteins and cholesterol-laden macrophages in the artery wall (Libby et al, 2019;Tabas et al, 2015).Certain macrophage behaviours are vital in enabling inflammation to resolve in both plaques and other inflammatory contexts.These behaviours include the emigration of macrophages from the site of inflammation and the clearance of dead cells (Tabas, 2010;Llodrá et al, 2004).Macrophages in advanced plaques, however, display an impaired ability to carry out normal anti-inflammatory functions (Tabas, 2010;Randolph, 2008).Understanding how these impaired behaviours may promote macrophage retention in early plaques is therefore important to identify the conditions that lead to advanced plaque development.In this paper, we develop a partial differential equation model for an atherosclerotic plaque based on a multiphase framework, and we use this to investigate how macrophage emigration and dead cell clearance influence the growth and composition of the early plaque.
Atherosclerotic lesions are initiated when cholesterol-carrying low density lipoproteins (LDL) from the bloodstream are deposited in the intima (Libby et al, 2019;Falk, 2006).The intima is an initially thin layer that separates the endothelium (a monolayer of endothelial cells lining the interior of the artery) from the media (a thicker layer comprising muscle cells and collagen).The endothelium becomes more permeable to particles such as LDL at sites of endothelial dysfunction, which occur due to disrupted blood shear flow or endothelial cell damage (Gimbrone Jr and García-Cardeña, 2013).LDL particles that have penetrated the endothelium and accumulated in the intima undergo oxidation and other forms of chemical modification to produce modified LDL (modLDL) (Madamanchi et al, 2005;Yoshida and Kisugi, 2010).Intimal modLDL is a potent immune trigger, and stimulates the expression of adhesion molecules by endothelial cells.These adhesion molecules bind to and capture monocytes circulating in the bloodstream, which then transmigrate through the endothelium into the intima (Blankenberg et al, 2003;Bobryshev, 2006).Once in the intima, monocytes differentiate into macrophages, which then take up modLDL via phagocytosis (Bobryshev, 2006;Moore et al, 2013).The resulting cholesterol-engorged macrophages are referred to as macrophage foam cells.Early atherosclerotic lesions consist largely of modLDL, macrophage foam cells, and debris from dead cells.
In atherosclerosis and other forms of inflammation, macrophages undergo a form of programmed cell death called apoptosis (Cohen and Mosser, 2013;Van Vré et al, 2012).Apoptosis is normal even in healthy tissue, and apoptotic cells express find-me and eat-me signals to encourage their clearance by live macrophages via efferocytosis.In chronic inflammatory conditions such as advanced atherosclerosis, the accumulation of apoptotic material will often overwhelm macrophages' efferocytic capacity.This is especially likely in cases where rates of macrophage death are elevated due to high levels of ingested cytotoxic material (such as cholesterol), or where efferocytosis itself becomes defective.Both of these are observed in atherosclerotic plaques (Yurdagul Jr et al, 2018;Schrijvers et al, 2005).Uncleared apoptotic cells will eventually undergo an uncontrolled form of cell death called necrosis.Necrotic macrophages are highly problematic due to their lower production of "find-me" and "eat-me" signals.Necrotic material is also highly inflammatory and will attract more macrophages which are themselves likely to undergo necrosis, thereby perpetuating the cycle.A major characteristic of advanced plaques is the presence of a large necrotic core, consisting of lipids and debris released by necrotic cells (Sakakura et al, 2013).The inefficient efferocytic clearance of apoptotic material is therefore an important precursor to the development of a vulnerable plaque state that is more likely to rupture.
The emigration of macrophages from sites of inflammation is a critical part of inflammation resolution in other inflammatory conditions (Lawrence and Gilroy, 2007;Bellingan et al, 1996), and is believed to play an important part in atherosclerosis as well, although its mechanisms are not yet well understood.Experimental studies have measured the expression of receptors such as CCR7 that are known to be involved in macrophage emigration in other inflammatory conditions, and find that higher expression of these receptors is correlated with reduced plaque size (Trogan et al, 2002;Finney et al, 2017).Several studies suggest that emigration most likely happens via migration into lymphatic vessels that connect to the outer artery.One such study found that after plaque regression was induced, macrophages that were originally present in the plaque were found in the lymph nodes (Llodrá et al, 2004).Another study involving fluorescent bead-tagged macrophages found that the distance between tagged macrophages and the elastic lamina (a membrane separating the intima from the outer muscle cell layers) decreased during plaque regression, suggesting that macrophages are exiting through the outer artery instead of transmigrating across the endothelium (Williams et al, 2018).
Mathematical modelling of the inflammatory response during atherosclerosis is an area of growing interest (Parton et al, 2016).Much of this work uses a continuum approach to investigate the activity of monocytes and macrophages in the vessel wall, and their interactions with lipoproteins.Previous work in our group has considered the accumulation and efferocytic clearance of apoptotic cells using both ODE models (Lui and Myerscough, 2021) and non-spatial lipid-structured PDE models for macrophage populations in plaques (Ford et al, 2019).An alternative model by Fok (2012) considers macrophages and dead cells in a reaction-diffusion model for plaque development.No other spatial models exist to date that explicitly consider both cell death and efferocytosis.Another shortcoming of many spatial models is their use of a fixed spatial domain to model the intima, where the domain remains fixed in size rather than expanding as the plaque develops.Some models employ a moving boundary to allow the artery wall to expand (Friedman and Hao, 2015;Thon et al, 2018b).
In particular, one model by Friedman and Hao (2015) employs a multiphase framework to capture the effects of cell crowding.However, none of these models incorporate macrophage emigration in a way that respects local mass conservation, although some fixed domain models incorporate an unbalanced sink term for macrophage removal due to death or emigration (Fok, 2012).
In this paper, we present a free boundary multiphase model for early plaque growth that focuses on the interplay between foam cell death, efferocytosis, and emigration.We formulate the model in Section 2. In Section 3, we investigate how the relative rates of foam cell death and efferocytic uptake influence plaque composition.In Section 4, we show how macrophage emigration can slow or halt plaque growth, and how its effectiveness depends on the availability of live foam cells in the deep plaque.In Section 5, we extend the model to include macrophage tagging, and we use the extended model to investigate macrophage retention in plaques.We conclude in Section 6 with a discussion of the results in the context of existing plaque models and biomedical studies, and some suggestions for future modelling work.

Model formulation
In this section we present a general model that describes the growth of an early stage atherosclerotic plaque.We model the intima as a 1-dimensional domain on the interval [0, R( t)], which represents a radial cross-section through the artery wall.The x = 0 boundary corresponds to the endothelium, and the x = R( t) boundary corresponds to the elastic lamina Figure 1: A 2-dimensional cross section of an artery, with a 1-dimensional radial cross section taken through the plaque for the problem domain.
that separates the intima from the media (see Figure 1).We allow the x = R( t) boundary to vary with time since early plaques cause the artery to undergo compensatory enlargement.During compensatory enlargement, the artery expands outwards in a way that preserves lumenal area (Korshunov et al, 2007).We also assume that the diameter of the artery is much larger than the width of the intima in early stage plaques (Bonithon-Kopp et al, 1996), so that Cartesian coordinates are a good approximation to radial coordinates.We notate all quantities with tildes for the full dimensional model, and then remove tildes for the nondimensional model in Section 2.1.
We assume that the plaque comprises three phases, each measured in cells or equivalent cell volumes per unit length.These phases are macrophage foam cells f (x, t), modified LDL (or modLDL) l(x, t), and dead cellular material c(x, t).For simplicity, we let f (x, t) encompass all monocyte-derived cells, including M1 and M2 macrophages and dendritic cells for instance (Murray and Wynn, 2011), ignoring differences in phenotypic expression.Similarly, we ignore the distinction between the various forms of chemically modified LDL, treating them as a single modLDL species l(x, t).For cell death, we ignore the distinction between apoptosis and necrosis, and use a generic cell death term to encompass both processes.
The phases obey the following continuity equations: In the above, Ju and su denote phase-specific flux and source terms respectively for the phases u = f, l, c, and ṽ is a common mixture velocity with which all phases are advected.
We also assume there are no voids throughout the plaque, so that the phases obey a no-voids condition f + l + c = N 0 (4) everywhere, where N 0 denotes the total phase density, which is constant with time and space.
Foam cells internalise modLDL via phagocytosis, and dead material via efferocytosis.We model both processes with a simple mass-action term.Foam cells also undergo cell death, which we model using a linear death term.This gives the following expressions for the source terms in equations ( 1) to (3): where μa is the macrophage foam cell death rate, μp is the rate of phagocytosis per foam cell per unit of available modLDL, and μc is the rate of efferocytosis per foam cell per unit of available dead material.Since sf + sl + sc = 0, the model is locally mass conservative.
We assume that foam cells undergo undirected random motion (Owen and Sherratt, 1997) and directed chemotactic motion towards modLDL.The latter is a simplification of the full process whereby modLDL stimulates the production of chemoattractants.Foam cells also undergo chemotactic motion towards dead material, in response to find-me signals expressed by the dying cells (Kojima et al, 2017).ModLDL and dead material undergo diffusive motion.For simplicity, we assume constant diffusion coefficients for all phases, and linear chemotaxis with regards to modLDL and dead material for foam cells.This gives the flux terms where Df is the foam cell random motility coefficient, Dl and Dc are diffusion coefficients for modLDL and dead material respectively, and χl and χc are chemotactic coefficients for foam cells in response to modLDL and dead material respectively.
The endothelium permits a net influx of LDL carried by native LDL particles, which we assume is constant over plaque growth timescales of days or weeks.We also assume that typical timescales for chemical modification of LDL are much faster than plaque growth timescales (Cobbold et al, 2002).Under this assumption, native LDL becomes inflammatory modLDL immediately after entering the intima, so we model native LDL deposition as an boundary modLDL influx.Monocyte recruitment requires adhesion molecules expressed by endothelial cells and chemokines produced in the intima.These are expressed in response to the presence of modified LDL.For simplicity, we assume monocyte recruitment is proportional to the amount of modLDL at the endothelium.We also assume that monocyte differentiation timescales are much faster than those for plaque growth (Blanchard et al, 1991), and so we model monocyte recruitment as a macrophage boundary influx.Cells do not die until they are already inside the intima, so the flux of dead material at the endothelial boundary is zero.This gives the following flux boundary conditions at x = 0: where σl is the rate of LDL deposition, and σf is the rate of monocyte recruitment per unit of modLDL at the endothelium.
At the medial boundary, we assume zero flux boundary conditions for the passive modLDL and dead material phases.We assume foam cells will emigrate into the lymphatic vasculature at a rate proportional to the quantity of foam cells present.The corresponding flux boundary conditions at x = R are given by where σ e is the average foam cell egress velocity.
For the initial conditions, we assume the intima is populated with a small number of live resident macrophages (Ley et al, 2011), so where x S is the diameter of a typical mouse macrophage.
In order to close the model, the mixture velocity ṽ(x, t) and the velocity of the medial boundary R ( t) need to be defined, in addition to initial conditions.Summing over the phase continuity equations ( 1) to (3) and applying the no-voids condition (4) gives Integrating equation ( 21) from x = 0 to x and substituting the x = 0 boundary conditions (11) to (13) gives the following expression for the mixture velocity: The boundary velocity can be similarly obtained by integrating equation ( 21) from x = 0 to R and substituting the x = 0 and R boundary conditions ( 11) to ( 16), which yields The full dimensional model consists of equations ( 1) to ( 20) and ( 23).

Parameter estimates and model nondimensionalisation
We nondimensionalise the model by applying the rescalings where t S and x S are some reference timescale and length scale, and N 0 is the total phase density.
We choose t S to be one week as experimental studies on plaques in mouse models typically take place over periods of weeks to months (Williams et al, 2018;Potteaux et al, 2011).
Lengths are rescaled by the diameter of a typical mouse macrophage.As observed macrophage sizes vary from under 14 µm for murine bone marrow-derived macrophages (Cannon and Swanson, 1992) to 20 µm for murine peritoneal macrophages (Nguyen et al, 2012), we choose a characteristic size of 16 µm.Estimating certain parameters requires the conversion of three-dimensional volume densities to one-dimensional linear densities.For the purposes of parameter estimation, we consider a cylindrical section through the intima whose circular cross-section is transverse to the positive x vector and has diameter x S , so that the maximum linear phase density N 0 is one cell per unit of x S .
Foam cell motility parameters Df and χl are estimated using models fitted to data from chemotaxis assay experiments (Owen and Sherratt, 1997;Sozzani et al, 1991).The phagocytosis parameter μp is based on experimental observations where macrophages are observed to phagocytose a significant volume of modLDL after a timescale of approximately 60 min (Zwaka et al, 2001).The boundary monocyte recruitment parameter σf is based on an experimental study by Jeng et al (1993) where monocytes pretreated with modLDL are exposed to shear flow and are observed to bind to an endothelial cell monolayer with an area of 0.1452 mm 2 at a rate of 140 -200 cells over 30 min.The death rate μa is based on experimental observations of mouse macrophages incubated in acetylated LDL, where 12% of the initial macrophage population is observed to undergo apoptosis after 9 h.Fitting this to a simple exponential decay curve gives an apoptosis parameter of about 0.25 h −1 before nondimensionalising.
The remaining parameters are based on order of magnitude estimates.The modLDL diffusion coefficient Dl is estimated to be larger than Df due to the small size of LDL particles compared to macrophages.Experimental studies of LDL phagocytosis (Zwaka et al, 2001;Suits et al, 1989) find that for heavily cholesterol-laden foam cells, cholesterol droplets comprise a significant proportion of the total foam cell volume, and so the range of LDL influx rates is chosen so that σl is of a similar order of magnitude to the monocyte influx σf l| x=0 .The upper range of foam cell egress velocities σe is chosen so that the medial egress flux σe f | x= R is of a similar order of magnitude to the endothelial LDL and monocyte fluxes σl and σf l| x=0 .For efferocytosis, we choose a range of values for μe for which efferocytic uptake rates are of a similar order of magnitude to apoptosis rates.Dead material diffusion is assumed to be slower than live macrophage random motility due to the lack of active locomotion.We choose Dc to be lower than Df within an order of magnitude.The macrophage chemotactic response to dead material χc is assumed to be lower than that for modLDL, since the dead material phase includes necrotic material which expresses few find-me signals.
Tables 1 and 2 summarise the model parameter values before and after rescaling.The final nondimensional model consists of the constitutive equations with no-voids condition source terms flux terms boundary conditions at x = 0 and initial conditions where the mixture velocity is and the domain grows at rate

Numerical solution
The system was solved numerically using the method of lines.The nondimensional system was first mapped onto a fixed spatial domain [0, 1] using the transformation (x, t) → (y, τ ) = (x/R(t), t), thereby transforming it into a system of parabolic PDEs coupled to an ODE for the domain length R(τ ).The transformed system was then reduced to a set of time-dependent ODEs by discretising spatial derivatives using a central differencing approximation.The resulting ODEs were solved using a backward differentiation formula method in Python using SciPy's integration libraries.
3 Plaque structure, cell death, and efferocytosis In this section we neglect emigration, setting σ e = 0, and focus on how the relative rates of cell death and efferocytosis influence the composition of the plaque.

Base case simulation
Figure 2 illustrates the evolution of a base case with relatively high efferocytic uptake compared to cell death.Following the initial injury where LDL starts being deposited in the intima, the plaque quickly reaches a state where it is populated almost entirely by foam cells and dead material, although the balance between the two may vary (Figure 2).A small amount of modLDL that has yet to be ingested remains near the endothelium, though some uningested modLDL may remain throughout the intima for very early times (Figure 2), or for cases where all foam cells in the deep plaque have died (Figure 3(a)).ModLDL and live foam cell densities are higher near the endothelium than in the deep plaque due to the constant influx of new monocytes and LDL particles, though some dead material is still present due to diffusion.Away from the endothelium, phase densities are largely uniform, with some non-uniformity observed near the medial boundary when σ e > 0 (Figure 3(d)). (a) x O D 9 K + f X 5 3 j c t y y a 6 V q p + q x Z P W 8 j T 4 7 H v j 6 P D N i v d e z i 4 8 9 J Z u 1 9 u l p Z R c S 3 o + Q F d J u t s q d 2 i 6 E r p B u t 9 6 1 r F 0 I X i F 2 t 9 5 o V X c h 7 g q p N S u d D 8 1 d y G S F V C p N q 9 x O D 9 K + f X 5 3 j c t y y a 6 V q p + q x Z P W 8 j T The intima itself quickly approaches a state where it grows at a constant rate (Figure 4).The rates at which foam cells and dead material accumulate also settle to a constant rate, while uningested modLDL may either reach a steady value or accumulate at a constant rate, depending on whether the deep plaque still has live foam cells remaining to ingest new lipid (Figure 4(a) vs (b)).
We remark that increasing the LDL deposition rate σ l or the monocyte recruitment parameter σ f will increase the intimal growth rate.This is due to the increased influx of new material in the form of modLDL or recruited monocytes respectively.As we see later however, the qualitative spatial structure depends primarily on the cell death and efferocytosis rates, which remain our focus in this section.
Figure 5(a) compares the relative sizes of the diffusive, chemotactic, and advective foam cell flux terms.Diffusion and chemotaxis are dominated by advective transport, which is the sole driver of foam cell movement in the deep plaque.Chemotaxis is important near the endothelium where there are larger amounts of uningested modLDL, but becomes irrelevant in the deeper plaque once the available modLDL has been internalised.From Figure 5(b), the total phase velocities for both the foam cell and dead cell phases converge to the mixture velocity in the deep plaque, with little interphase motion.This suggests that mass transport in the deep plaque occurs primarily via bulk advection.We note however that when σ e > 0, emigration of foam cells into the lymphatics causes some deviation between the mixture velocity and the foam cell phase velocity (Figure 5(c)).Nevertheless, the σ e = 0 case provides a good baseline for comparison, since the close agreement between the phase velocities makes the system amenable to analysis using the method of characteristics.

ODE approximation for deep plaque composition
The constitutive equations for each phase u can be written in the form where v u = (J u +vu)/u is the total velocity of phase u, and s u is its source term.In Figure 5, we observed that away from the endothelium, the bulk advection term vu dominates over random and chemotactic motion modelled by J u for all three phases, and so v u ≈ v in the deep plaque.For each phase u, we define U i (t) to be the phase density along a characteristic curve x c (t) advected with speed v u , i.e.
In the deep plaque, v u ≈ v for all phases, and we can make the approximation We note also that in the deep plaque, the modLDL density l(x, t) is negligible, so we make an additional approximation l = 0. Substituting the model source terms ( 29) and (30) into equation ( 49) yields the ODE system where F (t) = f (x c (t), t) and C(t) = c(x c (t), t) denote the phase densities along characteristic curves.Due to the no-voids condition (28), we also have F + C = 1.Using this condition, the coupled equations reduce to a single ODE This equation has two steady states at F * = 0 and F * = 1 − µa µe , with a transcritical bifurcation occurring when µ a = µ e .Linearising about the steady states gives eigenvalues of λ = µ e − µ a at F * = 0 and λ = µ a − µ e at F * = 1 − µa µe .When µ a < µ e and efferocytosis rates are high relative to cell death, F * = 1 − µa µe is the sole stable steady state, and the system tends towards a state where live foam cells coexist with dead material.When µ a < µ e however, F * = 0 is the sole stable state.Here, efferocytic recycling is insufficient to counterbalance cell death, and the plaque settles to a state where all foam cells die off.From the eigenvalues, for a fixed ratio µa µe , higher overall rates of cell death and efferocytosis will cause the system to settle to the relevant steady state faster.

Plaque composition and the death/efferocytosis balance
The deep plaque composition closely matches the steady state densities predicted by the ODE model.In Figure 3, lower rates of efferocytosis are associated with higher amounts of dead material and fewer living foam cells.For plaques with efficient efferocytosis where µ e > µ a (Figure 3(c) and (d)), the foam cell density settles to the f * = 1 − µa µe density predicted by the ODE model.On the other side of the µe µa = 1 bifurcation point, plaques with µ e < µ a (Figure 3(a) and (b)) reach a state where the deep plaque supports no live foam cells.For especially low values of µ e (Figure 3(a)), foam cells will die before ingesting all available extracellular modLDL, leaving some remnant modLDL which gets carried advectively into the deep plaque.From Figure 6, higher rates of death and efferocytosis for a given value of µa µe lead to better agreement between the ODE steady states and the deep plaque densities at x = R.This is consistent with linear analysis of the ODE system, which predicts faster convergence to the steady state for higher values of µ a and µ e .
Cell death and efferocytosis also influence the growth rate of a plaque, even with fixed values of the monocyte and LDL influx parameters σ f and σ l .From Figure 7(a), increasing µ e with fixed µ a will slow the rate of growth.This is because increasing the efficiency of efferocytosis will increase the numbers of foam cells actively consuming modLDL, thus reducing endothelial modLDL levels and dulling the inflammatory response.This slowing in growth can also be observed when comparing plaque sizes in Figure 3. Reducing µ e and µ a with constant µe µa has the same effect of slowing plaque growth by reducing endothelial modLDL levels, since reducing the rate at which cells start to die will increase the quantity of available foam cells near the endothelium.
The rate at which a plaque accumulates new dead material is a useful marker of plaque health alongside its total growth rate.Both are low in a healthy plaque.As discussed previously, reducing µ e and µ a with constant µe µa reduces the rate at which new cellular material enters the plaque.From Figure 7(b), this results in a slight decrease in the total rate of increase of dead material, although the proportion of dead material in the bulk of the plaque remains the same.Increasing the efficiency of efferocytosis via µ e with constant µ a has a more significant effect on the rate of accumulation, as it reduces both the monocyte influx and the local dead cell proportion in the interior of the plaque by promoting more active uptake of dead material.

Emigration and plaque growth
Foam cell emigration can slow plaque growth and potentially allow plaque size to stabilise by allowing material to leave the plaque.The ability of emigration to impact plaque growth however depends on the availability of live foam cells.
In Figure 3(a,b) and 8(a), µ e < µ a , and efferocytosis is insufficient to prevent the death of foam cells in the deep plaque.Increasing the emigration velocity here has little effect on plaque size as there are few foam cells close to the medial edge of the plaque.In Figure 3(c,d) on the other hand, µ e > µ a , and there is a higher density of foam cells in the deep plaque.Increasing the emigration velocity here slows the rate of plaque growth by enabling material to leave the plaque, which also results a reduction in plaque size.For sufficiently high σ e , the plaque will stabilise entirely and settle to a steady state with no growth (Figure 8(b)).While emigration does cause a slight reduction in foam cells numbers near the media, their numbers remain sufficient to allow emigration.We remark briefly that emigration causes phase densities near x = R to deviate from those predicted by the ODE model, which we discuss further in Section 6.
The importance of healthy efferocytosis in enabling foam cell emigration is more apparent when considering a range of values of µ e .From Figure 9(b), increasing the emigration velocity σ e will increase the total foam cell emigration flux for most values of µ e .This flux is diminished for lower values of µ e , and for µ e ≤ µ a , foam cell emigration remains negligible even for very high values of σ e .Figure 9(c) suggests that for cases where when µ e ≤ µ a , the lack of emigration observed is due to the absence of live foam cells in the deep plaque.Conversely, cases with µ e > µ a will still have some foam cells near the media, thus enabling emigration.Even for these cases however, higher emigration velocities result in lower medial foam cell densities due to the faster removal of foam cells.In particular, for cases where µ e > µ a with µ e close to µ a , the low live foam cell density gives diminishing returns on the total emigration flux when the emigration velocity is increased.
These results suggest that cell death will interfere with or prevent plaque stabilisation via foam cell emigration, depending on how efficiently efferocytosis acts.From Figure 9(a), for µ e ≤ µ a , increasing σ e has no effect on the plaque growth rate as there are no foam cells to emigrate from the plaque.For values of µ e > µ a close to µ a , emigration is able to slow plaque growth, but full plaque stabilisation is not observed even for very high values of σ e .For higher values of µ e , stabilisation is observed for high enough σ e , with higher µ e scenarios achieving this more readily.We briefly note that the endothelial LDL and monocyte influx parameters σ l and σ f also affect plaque growth in the presence of emigration.This is because higher σ l and σ f necessitate higher emigration rates σ e to balance the higher rates of material influx at the endothelium.

Macrophage circulation and retention
In this section, we investigate the movement and retention of macrophage foam cells in the intima.We do so by extending the model to include an additional bead species that moves with foam cells (or dead cells), based on experimental work that labels and tracks macrophages using fluorescent latex microspheres (Williams et al, 2018).

Extending the model: bead tagging
For the bead labelling model, we add two new species q f and q c , representing beads carried by live foam cells and dead material respectively.We assume that the volume occupied by the beads is negligible, and that their densities can be modelled by the continuity equations We suppose that beads are transported with the same velocity as the phase in which they are being carried.Using the flux-velocity relation v u u = j u (= J u + vu) gives the bead flux terms Beads carried by macrophages are released into the dead material phase at the rate µ a with which foam cells die.Similarly, beads within dead material are taken up by foam cells via efferocytosis at the same rate µ e m as dead material is taken up.This gives the bead source terms Note that bead numbers are locally conserved since s q f + s qc = 0.
At the endothelial boundary, we assume that macrophage-carried beads are introduced over a brief time period centred at t tag .We model the bead influx with a Gaussian function, so that bead counts rise and fall over approximately 6 hours (t = 0.035).We represent this via the following bead endothelial boundary conditions: The corresponding medial boundary conditions are These equations can be derived from equations ( 16) and (55).

Results
For plaques with sufficiently high foam cell emigration, foam cells will move closer to the medial boundary and eventually reach it, before emigrating into the arterial lymphatics.
Figure 10 and Figure 11(a) show the time evolution of the bead distribution in a plaque with nonzero emigration (σ e = 20) and efferocytosis rates that allow live foam cells to persist deep into the plaque (µ a = 40, µ e = 80, so µe µa = 2).After bead-carrying monocytes have been recruited, the bead population is pushed deeper into the plaque as the plaque grows.Despite the increasing plaque size, the distance from the bead population to the medial boundary decreases due to material exiting through the boundary.
The amount of time that foam cells spend in the intima depends on emigration rates and the size of the plaque.Figure 11 considers the same high efferocytosis plaque scenario as Figure 10, but with bead tagging happening at different times.From Figure 11(a), as the plaque increases in size, tagged foam cells that are recruited later have to travel a greater distance before they reach the medial boundary.As a result, they spend more time in the plaque.Figure 11(b) shows the time evolution of the total amount of beads Bead totals peak soon after the tagging time t tag , and fall to 0 when tagged foam cells emigrate from the intima.Bead-tagged foam cells that are recruited later spend more time in the intima.Figure 12 shows how the total circulation time of tagged monocytes depends on the time of recruitment and the average emigration velocity.The total circulation time t circ is estimated using the temporal full width at half maximum duration of the peak bead quantity, i.e.
Tagged monocytes that are recruited later circulate for longer time periods.Increased migratory behaviour of foam cells leads to reduced circulation times.This is because higher values of σ e slow plaque growth, reducing the plaque size and the distance foam cells have to travel before exiting the plaque.As observed earlier in Figure 5(b) and (c), higher values of σ e also increase the velocity of foam cells near the medial boundary.
For plaques with poor efferocytic uptake rates that contain mostly dead cells, bead-tagged foam cells will die, and the beads will remain in the plaque indefinitely.Figure 13(a) tracks the evolution of the bead distribution for a plaque with nonzero emigration velocity (σ e = 40) and poor efferocytosis (µ a = 40, µ e = 80, so µe µa = 0.5).Like the high efferocytosis case, bead-carrying cells are pushed deeper into the plaque as new material enters through the endothelium.Here, however, beads remain equidistant from the medial boundary as the plaque grows.Figure 13(a) tracks bead totals Q(t) with time, and shows that bead numbers peak soon after the monocyte tagging time t tag , and continue to remain high at longer times.This suggests that bead-carrying cells are failing to clear from the plaque, as no live foam cells are available to emigrate and allow material to leave the intima.

Discussion
In this work, we have constructed a new model for early atherosclerotic plaque growth that includes the effects of cell death on plaque development.A key aspect of our approach is the  ).Plot (a) shows the time evolution of the total bead distribution q f (x, t) + q c (x, t), where beads are administered multiple times at t tag = 0.5, 1.0, 1.5, 2.0, and 2.5.Plot (b) compares the time evolution of the total amount of beads Q(t) = R(t) 0 (q f + q c ) dx for each separate bead dosage.use of a locally mass-conservative multiphase framework with a free boundary.This allows us to model growth of the intima due to the accumulation of new material, and also to model cell movement and emigration in a manner that respects local mass conservation.
Our model finds that dead cells are most plentiful in the deep plaque, with a higher proportion of live macrophages closer to the endothelium.This is consistent with histological studies of early human aortic fatty streaks (Guyton and Klemp, 1993), which contain macrophage-derived foam cells mostly in their middle and closer to the endothelium.Deposits of cholesterol crystals are observed deeper inside the intima, which are indicative of foam cell necrosis.Observations of more advanced plaques in humans and mice show a similar structure, where necrotic cores form deep in the plaque, and live macrophages tend to localise closer to the endothelium, albeit with a much larger necrotic core than in early fatty streaks (MacNeill et al, 2004;Dickhout et al, 2011).
Our model also suggests that advection is the primary mass transport mechanism that pushes material deeper into the plaque.The multiphase framework we use is highly idealised, where we include a simple mixture velocity advection term to help enforce the no-voids condition.This is in contrast to a force-balance equation approach or a full continuum mechanical approach based on porous flow theory or fluid mechanics, where advection-type terms will arise indirectly from other assumptions.Nevertheless, the mixture velocity advection term here plays a similar role to the advection term that arises in porous flow plaque models (Friedman, 2015;Cilla et al, 2015), which is typically induced by transmural pressure gradients across the endothelium or intima .
Despite the necrotic core being a key feature of advanced plaques, the role of dead cells in plaque development is largely neglected by other spatial models, and those spatial models that include cell death tend to model it via a simple unbalanced sink term.Ours is also the first spatially structured plaque model to explicitly consider efferocytosis, although Ford Figure 13: Bead circulation for different bead tagging times for a plaque with poor efferocytosis (µ a = 40, µ e = 30, σ e = 40).Plot (a) shows the time evolution of the total bead distribution q f (x, t) + q c (x, t), where beads are administered multiple times at t tag = 0.5, 1.5, and 2.5.Plot (b) compares the time evolution of the total bead quantity Q(t) = R(t) 0 (q f + q c ) dx for the three bead dosages in (a).et al (2019) consider efferocytosis in a lipid-structured model without spatial structure.Fok (2012) presents a reaction-diffusion-chemotaxis model for necrotic core formation that explicitly includes dead cells.His model posits that the accumulation of macrophages deep in the plaque is driven mostly by chemotaxis towards extracellular lipid in the core, whereas in our model, advection is the primary means of mass transport into the deep plaque.Fok's model also neglects the recycling of efferocytosed material by live cells, although it includes a term for the removal of dead cells by other efferocytosis-competent leukocytes.
In our model, the rate of efferocytic recycling relative to cell death is a key determinant of whether the deep plaque can support live foam cells.Higher rates of efferocytosis relative to death result in a higher proportion of live cells.These results are consistent with gene knockout studies in mice, which find that defective efferocytosis is correlated with large plaque size and higher numbers of dead cells (Kojima et al, 2017).Crucially, our model undergoes a bifurcation as the ratio of the efferocytosis and death parameters changes, and complete foam cell death in the deep plaque is unavoidable if the ratio falls below a threshold value.Studies have explored the potential for using efferocytosis as a therapeutic target to treat cardiovascular disease and cancer (Kojima et al, 2017), with strategies including disruption of the expression of "don't-eat-me" signalling molecules (Kojima et al, 2019), and the use of pro-efferocytic nanoparticles (Flores et al, 2020) for example.Our model has the unfortunate implication that small improvements in macrophage efferocytic uptake may not be enough to induce regression of unstable advanced plaques.Larger improvements however may enable regression by increasing live foam cell counts.
Foam cell emigration in our model is critical in slowing plaque growth, and can allow plaques to stabilise if emigration is large enough.This is consistent with experimental studies in mice in which the CCR7 receptor is blocked.CCR7 is involved in macrophage emigration, and blocking the receptor leads to larger plaques with impaired macrophage migration (Luchtefeld et al, 2010).Intuitively, model plaques that are larger and plaques with slower foam cell emigration are characterised by increased retention of foam cells.This has implications for necrotic core formation, since foam cells that emigrate sooner are less likely to become necrotic.Emigration has been observed to enable steady states in other models, including the lipid-structured model by Ford et al (2019), and in ODE models by Cohen et al (2014) and Lui and Myerscough (2021).Emigration is typically neglected in existing spatial models however, and those that do exhibit steady states typically model emigration using unbalanced macrophage sink terms, rather than via emigration through the domain boundaries.
Inefficient efferocytosis and high rates of cell death inhibit foam cell emigration and plaque stabilisation by lowering the numbers of live foam cells.Below the threshold efferocytosis/death ratio required for the survival of deep plaque foam cells, emigration is effectively negligible.High rates of macrophage death and poor rates of emigration are both observed in plaques in hypercholesterolemic mice (Feig et al, 2011).However, it remains unclear whether the lack of live cells specifically contributes to poor emigration rates in vivo, or whether both are a consequence of other factors such as excessive cholesterol loading.Attempts to block apoptotic pathways (using caspase inhibitors for instance) have been inconclusive due to the subsequent increase in necrosis in cells with dysfunctional apoptotic pathways (Martinet et al, 2011).
Our model provides a good foundation for developing more detailed models of early plaque growth.One immediately relevant question is the effect of intracellular cholesterol, which is known to upregulate apoptotic pathways and disrupt normal cell function in macrophages (Tabas, 2002).By separating the foam cell phase into two individual phases representing macrophages and intracellular cholesterol, we can model heterogeneity in foam cell cholesterol loads.The intracellular cholesterol phase would inherit the macrophages' phase velocity, like the bead phase in this paper.This would allow us to incorporate some of the advantages of lipid-structured PDE models (Meunier and Muller, 2019;Ford et al, 2019;Watson et al, 2022) into a spatial framework, such as cholesterol-dependent death and emigration rates.
Another important extension is the inclusion of HDL, which has multiple atheroprotective functions (Barter et al, 2004;Brites et al, 2017;Feig et al, 2011), including the ability to accept cholesterol from foam cells, inhibit LDL oxidation, and stimulate macrophage emigration.Other modelling work has shown that increasing HDL levels can enable plaque stabilisation (Thon et al, 2018a;Cohen et al, 2014;Ford et al, 2019;Friedman and Hao, 2015).By extending our model to include HDL, we can explicitly consider how it reduces the accumulation of apoptotic and necrotic material.
The inclusion of a bead cargo species in our model also provides a promising means to explore further questions raised by Williams et al (2018) on macrophage motility and retention in plaques.Tracking individual cells is trivial in discrete cell-based models, but presents a challenge in a continuum modelling framework, and our bead model represents a novel method for tracking the movement of cells within a spatial PDE model.This approach may also prove useful in modelling other diseases where macrophage motility is an important question, such as in granulomas (Egen et al, 2008) and tumour infiltration (Owen et al, 2004).

Figure 5 :
Figure 5: (a) Comparison of the individual flux terms in the constitutive equation (equations (25) and (29))) for the foam cell phase: random motion −D f ∂f ∂x , chemotaxis towards modLDL χ l f ∂l ∂x and dead material χ c f ∂c ∂x , and bulk advection vf .(b,c) Comparison of the foam cell and dead material phase velocities with the mixture velocity.Terms are plotted as a function of x at a fixed time t = 1 for a scenario with σ a = 40 and (b) σ e = 0, (c) σ e = 80.

Figure 6 :
Figure 6: Comparison of the deep plaque foam cell density f | x=R as a function of the efferocytosis/death ratio µe µa , plotted at t = 0.5 for various values of µ a .The equilibrium density predicted by the ODE model is also given for comparison.

Figure 7 :
Figure 7: Comparison of (a) the intima growth rate dR dt and (b) the rate of change of total dead material d dt

Figure 9 :
Figure 9: Plots showing the effect of the emigration velocity σ e on plaque development, compared across various efferocytosis rates µ e with µ a = 40 fixed.Quantities are measured at t = 1.Plots compare (a) the total intima growth rate dR dt , (b) the total flux of foam cells out of the medial boundary j f | x=R via emigration, and (c) the foam cell density at the medial boundary f | x=R .

Figure 10 :
Figure 10: Time evolution of the bead spatial distribution for a plaque with high efferocytosis (µ a = 40, µ e = 80, µ e = 20), where the bead tag dosage is centred at time t tag = 1.0.Distributions are plotted at times (a) t = 1.1,(b) t = 1.5, and (c) t = 1.8.The solid black line denotes the medial plaque boundary x = R(t).

Figure 11 :
Figure11: Bead circulation for different bead tagging times for a plaque with high efferocytosis (µ a = 40, µ e = 80, σ e = 20).Plot (a) shows the time evolution of the total bead distribution q f (x, t) + q c (x, t), where beads are administered multiple times at t tag = 0.5, 1.0, 1.5, 2.0, and 2.5.Plot (b) compares the time evolution of the total amount of beads Q(t) =

Figure 12 :
Figure 12: Total bead circulation time t circ for varying bead tagging times t tag , plotted for different emigration velocities σ e .Plaque parameters used are µ a = 40, µ e = 80.

Table 1 :
Dimensional parameter estimates (order of magnitude estimates used where quantitative biological data was unavailable)

Table 2 :
Rescaled parameters for the nondimensionalised model