A Lipid-Structured Model of Atherosclerotic Plaque Macrophages with Lipid-Dependent Kinetics

Atherosclerotic plaques are fatty growths in artery walls that cause heart attacks and strokes. Plaque formation is driven by macrophages that are recruited to the artery wall. These cells consume and remove blood-derived lipids, such as modified low-density lipoprotein. Ineffective lipid removal, due to macrophage death and other factors, leads to the accumulation of lipid-loaded macrophages and formation of a necrotic lipid core. Experimental observations suggest that macrophage functionality varies with the extent of lipid loading. However, little is known about the influence of macrophage lipid loads on plaque fate. Extending work by Ford et al. (J Theor Biol 479:48–63, 2019) and Chambers et al. (A lipid-structured model of atherosclerosis with macrophage proliferation, 2022), we develop a plaque model where macrophages are structured by their ingested lipid load and behave in a lipid-dependent manner. The model considers several macrophage behaviours, including recruitment to and emigration from the artery wall; proliferation and apotosis; ingestion of plaque lipids; and secondary necrosis of apoptotic cells. We consider apoptosis, emigration and proliferation to be lipid-dependent and we model these effects using experimentally informed functions of the internalised lipid load. Our results demonstrate that lipid-dependent macrophage behaviour can substantially alter plaque fate by changing both the total quantity of lipid in the plaque and the distribution of lipid between the live cells, dead cells and necrotic core. The consequences of macrophage lipid-dependence are often unpredictable because lipid-dependent effects introduce subtle, nonlinear interactions between the modelled cell behaviours. These observations highlight the importance of mathematical modelling in unravelling the complexities of macrophage lipid accumulation during atherosclerotic plaque formation.


Introduction
Atherosclerotic plaques are localised accumulations of cells, lipids and associated debris that form in the lining of major arteries (Hansson and Libby 2006). Plaques are initiated when blood-borne low-density lipoprotein (LDL) penetrates the endothelium of the artery and is deposited in the artery wall (Tabas et al. 2007). Once inside the artery wall, LDL is oxidised or modified in different ways. Accumulation of modified LDL elicits an immune response that attracts circulating monocytes. These monocytes rapidly differentiate into macrophages, which ingest (phagocytose) the modified LDL and stimulate further macrophage recruitment through inflammatory signalling (Moore et al. 2013;Tall and Yvan-Charvet 2015). The death of lipid-loaded macrophages (known as foam cells) creates fatty deposits that may accumulate over time to form a large necrotic core (Lusis 2000). Rupture of a vulnerable plaque can release this necrotic material into the bloodstream and trigger a clotting cascade that causes stroke or myocardial infarction (Lusis 2000;Hansson and Libby 2006).
Not all plaques progress to become clinically dangerous. Many simply resolve naturally or evolve towards a benign, non-resolving state (Bäck et al. 2019). The fate of a plaque is largely determined by the interaction between macrophages and lipids in the artery wall and, in particular, the relative rates at which these constituents enter and leave the tissue (Moore et al. 2013). In addition to recruitment by inflammatory signalling, the plaque macrophage population can be increased by local proliferation (Robbins et al. 2013;Lhoták et al. 2016). On the other hand, plaque macrophage numbers can be reduced by death (apoptosis) or by emigration out of the wall (Tabas 2010;Llodrá et al. 2004). Since infiltrating LDL becomes bound to the artery wall extracellular matrix (Tabas et al. 2007), the removal of this lipid from the system requires the intervention of macrophages. Lipid internalised by macrophages can be ferried out of the plaque during macrophage emigration or by macrophage offloading to high-density lipoprotein (HDL) in a process known as reverse cholesterol transport (Yvan-Charvet et al. 2010). A further important mechanism, which recycles the cell and lipid content of the plaque, is macrophage efferocytosis (Yin and Heit 2021). As observed in vitro, a macrophage can phagocytose an entire apoptotic cell and acquire the dying cell's ingested lipid (Ford et al. 2019b). Apoptotic cells that are not efficiently cleared by efferocytosis are a source of necrotic material (Kojima et al. 2017).
Experimental results indicate that plaque macrophage behaviour can change with the extent of lipid loading (Tabas and Bornfeldt 2016). Lipid accumulation in plaque macrophages can upregulate the production of the pro-inflammatory signals that are required for monocyte recruitment (Tall and Yvan-Charvet 2015). The cytotoxic effects of excessive lipid ingestion can lead to macrophage apoptosis (Tabas 2002;Feng et al. 2003). Using a detailed transcriptomic study of murine plaque macrophages, Kim et al. (2018) showed that macrophage proliferation is likely to decrease with increasing lipid load. Results from both in vitro and in vivo studies further suggest that macrophages with larger lipid loads are less likely to emigrate from plaques (Chen et al. 2019;van Gils et al. 2012;Wanschel et al. 2013). This is due to either reduced migration capacity (Chen et al. 2019) or through increased expression of so-called retention factors (van Gils et al. 2012;Wanschel et al. 2013). Defective macrophage efferocytosis, an often cited mechanism of atherosclerotic plaque progression (Thorp and Tabas 2009;Linton et al. 2016), may also arise through lipid-dependent effects; lipid loading can reduce the efficiency of efferocytosis by disrupting relevant signalling pathways (Yin et al. 2020) and excessive lipid acquisition through efferocytosis can lead to cytotoxic macrophage death (Yin and Heit 2021).
Experimental plaque formation data is frequently collected ex vivo (e.g., from studies that use the genetically modified, atherosclerosis prone ApoE −/− mouse). Thus, although the lipid-dependent effects noted above have been identified, little is known about how they contribute either individually or collectively to plaque formation dynamics. In this paper, we address this issue by developing a mathematical model to study the impact of lipid-dependent macrophage behaviour on the dynamics and fate of atherosclerotic plaque progression. As will be seen, our findings indicate that lipid-dependent effects can substantially alter the fate of a plaque, often in unpredictable ways. This provides valuable new insight, and emphasises the need to consider lipid-dependent plaque cell behaviour both in future mathematical models of plaque formation and in the interpretation of experimental plaque formation data.
Interest in mathematical modelling of atherosclerotic plaque formation has grown in recent years (Parton et al. 2016;Avgerinos and Neofytou 2019;McAuley 2021). Most work to date has focussed on modelling the inflammatory response of macrophages in the early plaque. Published approaches include spatially-averaged ODE models (Bulelzai and Dubbeldam 2012;Cohen et al. 2014;Islam and Johnston 2016;Thon et al. 2018;Lui and Myerscough 2021), spatially-resolved PDE models (El Khatib et al. 2007;Calvez et al. 2009;Little et al. 2009;Yang et al. 2016;Chalmers et al. 2017; Thon et al. 2019;Silva et al. 2020) and agent-based models (Bhui and Hayenga 2017). It is common to assume that the modelled macrophages have two sub-populations: those with little or no internalised lipid (usually termed macrophages) and those with lots of internalised lipid (usually termed foam cells) (Calvez et al. 2009;Bulelzai and Dubbeldam 2012;Cilla et al. 2014;Hao and Friedman 2014;Islam and Johnston 2016;Yang et al. 2016;Chalmers et al. 2017;Silva et al. 2020). Lipid-dependent macrophage behaviour can be implicitly incorporated in this model framework by assuming that these sub-populations have, for example, different rates of lipid consumption (Calvez et al. 2009;Bulelzai and Dubbeldam 2012;Cilla et al. 2014;Hao and Friedman 2014;Islam and Johnston 2016;Silva et al. 2020), migration (Calvez et al. 2009;Cilla et al. 2014;Yang et al. 2016;Chalmers et al. 2017;Silva et al. 2020) or apoptosis (Hao and Friedman 2014;Islam and Johnston 2016;Silva et al. 2020). An alternative approach is to model macrophages as a single population and track the population's total ingested lipid content (Little et al. 2009;Cohen et al. 2014; Thon et al. 2018Thon et al. , 2019Lui and Myerscough 2021). Here, lipid-dependent effects can be included at the populationlevel by assuming that macrophage behaviour depends on the average ingested lipid load (Thon et al. 2018(Thon et al. , 2019. A more natural means to model lipid accumulation in plaque macrophage populations, including with lipid-dependent effects, is to use a population model in which macrophages are physiologically structured by their internalised lipid content. Lipid-structured models of plaque macrophages have been developed by Ford et al. (2019a), Chambers et al. (2022), and Meunier and Muller (2019).
The model by Ford et al. (2019a) uses a system of partial integro-differential equations to study how internalised lipid loads are distributed in live and apoptotic plaque macrophage populations, and how this influences necrotic core formation. Possible behaviours of live macrophages in the model include: (i) apoptosis; (ii) emigration from the plaque; (iii) LDL and necrotic lipid phagocytosis; (iv) lipid offloading to HDL; and (v) efferocytosis of apoptotic cells. Simulations and analysis of this model demonstrate the important roles of emigration and efferocytosis in the prevention of necrotic core growth. Efferocytosis is identified as a double-edged sword, however, as it can drive ingested lipid loads to become excessively large [as confirmed experimentally in Ford et al. (2019b)]. Chambers et al. (2022) extended the Ford model to include macrophage proliferation. This introduces an additional means of reducing cell lipid loads as internalised lipid in the parent cell is split between its daughter cells upon division. The model demonstrates that macrophage proliferation can reduce necrotic core formation by enhancing the capacity for necrotic lipid consumption. However, results suggest that proliferation can also be a double-edged sword because the reduction in necrotic core often comes at the expense of a substantially enlarged macrophage population. Ford et al. (2019a) and Chambers et al. (2022) assume that all plaque macrophage behaviours occur at constant rates independent of internalised lipid content. In this paper, we generalise their lipid-structured framework to include live macrophage behaviour that depends smoothly and continuously on ingested lipid load. For now, we include lipid-dependent behaviour only in macrophage apoptosis, emigration and proliferation. Modelling of lipid-dependent phagocytosis and efferocytosis will be addressed in a future study. By simulating and analysing this new model, we demonstrate that lipid-dependent macrophage behaviour can substantially alter plaque fate by changing both the distribution and the net accumulation of lipid in the system. Note that, while lipid loading is believed to influence plaque fate by modulating the phenotypic profile of the macrophage population (Moore et al. 2013;Tabas and Bornfeldt 2016), we do not explicitly consider macrophage phenotype here. The current work can be regarded as a step towards more detailed mathematical models that link macrophage phenotype to internalised lipid load.
The remainder of the manuscript is structured as follows. In Sect. 2, we outline our methodology. This includes a brief presentation of the model equations, and definitions of the functions that characterise lipid-dependent macrophage behaviour. Section 3 reports results and analysis from an in-depth simulation study that addresses how lipid-dependent apoptosis, emigration and proliferation influence plaque progression. We discuss the implications of our results for both theoretical and experimental atherosclerosis research in Sect. 4, and end with broad conclusions on the significance of the study in Sect. 5.

Definitions
The number densities of live and apoptotic plaque macrophages with lipid load a a 0 at time t 0 are denoted m (a, t) and p (a, t), respectively. The minimum lipid load a 0 represents the endogenous lipid contained in the internal structures of each cell. We denote the acellular necrotic lipid content of the plaque at time t by N (t). For m (a, t) and p (a, t), we define the total number of cells in each population by: and the total lipid content of each population by: respectively. Lipid-dependent cell behaviour is modelled by assuming that the dimensional reference rate of a given behaviour is modulated by a dimensionless factor g (a). The symbol is a placeholder that represents β for apoptosis, γ for emigration or ρ for proliferation. For notational convenience, we define the related quantities: which will appear in the model equations below. Lipid-independence in a given behaviour can be recovered by setting g (a) ≡ 1. This, in turn, gives G (t) = M (t) and G a (t) = A M (t).

Model Statement
A formal derivation of the original lipid-independent model can be found in Ford et al. (2019a) and, for the macrophage proliferation terms, in Chambers et al. (2022). The interested reader is directed to these works for in-depth explanation. Here, we give a brief summary of the modelling assumptions before stating the equations in full. The model considers the following processes in the plaque: 1. Live macrophages consume lipid from LDL particles and offload lipid to HDL particles. Lipid attached to LDL particles is assumed to enter the plaque at constant rate λ L and be rapidly consumed by macrophages. A quasi-steady state approximation gives λ L M as the rate per cell of lipid ingestion from LDL. HDL particles are assumed to enter the plaque with a constant lipid capacity per time λ H and rapidly fill with lipid offloaded from macrophages. A quasi-steady state approximation gives λ H M as the rate per cell of lipid offloading to HDL. The net rate of lipid acquisition per cell from these two processes is λ L −λ H M = λ M . We assume that λ > 0 is necessary for a plaque to form. This continuous acquisition of lipid is modelled by an advection term in the m (a, t) equation. 2. Live macrophages consume necrotic lipid at a rate proportional to θ . This phagocytic lipid uptake is also modelled by continuous advection in the m (a, t) equation. 3. Live macrophages become apoptotic macrophages at rate βg β (a) or emigrate from the plaque at rate γ g γ (a). 4. Live macrophages consume apoptotic macrophages by ingesting their entire lipid content at a rate proportional to η. This process (efferocytosis) contributes a local sink term and a non-local source term to the m (a, t) equation. The source term takes the form of a convolution integral, which accounts for all possible consumption events that produce live macrophages with lipid content a. The value of this integral is interpreted as 0 for all a 2a 0 . 5. Live macrophages proliferate at rate ρg ρ (a). Shortly before division, we assume that the parent cell newly-synthesises the quantity of endogenous lipid a 0 that is required to form a second daughter cell (Scaglia et al. 2014). Upon division, we assume that the total lipid content of the parent is divided equally between the two daughters. Proliferation contributes a local sink term and a non-local source term to the m (a, t) equation. The source term accounts for the fact that daughter cells with lipid content a are produced by parent cells with lipid content 2a − a 0 . The transient, pre-division increase in the parent cell lipid content is not explicitly modelled as it occurs on a timescale much shorter than the timescale of interest. 6. Apoptotic macrophages undergo post-apoptotic necrosis, producing necrotic lipid at rate ν. 7. Live macrophages with lipid content a 0 are recruited to the plaque from the bloodstream. The rate of cell recruitment is assumed to be a saturating function of the total ingested lipid in the live macrophage population A M (t) − a 0 M (t), with maximal recruitment rate α and half-maximal recruitment when A M (t)−a 0 M (t) = κ. Note that this model for macrophage recruitment encodes an implicit lipid-dependence because the expression A M − a 0 M arises from an assumption that macrophages produce recruitment-stimulating cytokines at a rate proportional to their accumulated lipid content a − a 0 . Recruitment is modelled by a boundary condition on the m (a, t) equation.
The equations and boundary condition that reflect the above assumptions are as follows: λ Additional ODEs for M (t), A M (t), P (t) and A P (t) can be generated by integrating (4) and (5) with respect to a, both with and without multiplying by a prior to the integration. For these calculations, we impose the requirement that m (a, t), am (a, t), p (a, t) and ap (a, t) all → 0 as a → ∞. This leads to the following equations: Based on a similar argument to that presented in Chambers et al. (2022), we note that Eqs. (8) and (9) will lead to unbounded growth of M (t) and A M (t) if: To avoid this eventuality in the current study, we ensure that the condition: is satisfied for all a ∈ [a 0 , ∞).
In the presence of lipid-dependent cell behaviour, the ODE equations (6), (8)-(11) and the PDE equations (4)-(5) cannot be readily decoupled, in contrast to the simpler models of Ford et al. (2019a) and Chambers et al. (2022). This limits the opportunity for analytical investigation of the model equations and, hence, this paper will use numerical simulations to study plaque fate and dynamics in a range of scenarios.
The model is closed by assigning appropriate initial conditions to each variable. The PDE variables require initial distributions and the ODE variables require initial values. Generically, we set: For m 0 (a) and p 0 (a), we use the following half-normal distributions: which are scaled such that ∞ a 0 m 0 (a) da = M 0 and ∞ a 0 p 0 (a) da = P 0 . The parameter a σ > 0 defines the shape of the distributions. For A M0 and A P0 , we correspondingly define: The initial number of live macrophages M 0 can be defined in terms of a σ by assuming that the initial distribution m 0 (a) satisfies the boundary condition (7) at t = 0. This leads to the relationship: For a valid (positive) M 0 value, we require a σ > λ √ 2 α √ π . Finally, we assume that initially there is no necrotic lipid in the system (N 0 = 0), and the system contains fewer dead cells than live cells. To satisfy the latter assumption, we arbitrarily set P 0 = 0.5M 0 .

Lipid-Dependent Rate Functions
In the model, lipid-dependent cell behaviour is described by either a monotonic function g (a) = g s (a) or a non-monotonic function g (a) = g r (a). For the monotonic function, we use the following saturating relationship: Here, g s (a 0 ) = 1 and the non-negative parameter δ = lim a→∞ g s (a). Hence, g s (a) increases from 1 to δ when δ > 1 or decreases from 1 to δ when 0 δ < 1. The exponent n 1 determines the shape of the function and the parameter a > a 0 denotes the lipid content for which g s (a ) = 1+δ 2 . For the non-monotonic function, we use the following scaled relationship: where 0 < 1, b > 0 and the exponents q and k satisfy q > k 1. This function increases from g r (a 0 ) = to a peak value of 1 at a = a 0 + b q k q −k , before decreasing towards as a tends to infinity. The exponent k controls the rate of increase of g r (a) for a ≈ a 0 , while the larger exponent q controls the rate of decline of g r (a) as a → ∞.
In practice, we shall primarily focus on modulating factors that have a monotonic dependence on a (i.e. monotonic increasing for g β (a), monotonic decreasing for g γ (a) and g ρ (a)). However, for g γ (a), we shall also consider the non-monotonic dependence, since physical arguments suggest that this may be a more realistic assumption.

Nondimensionalisation
Using tildes to denote dimensionless quantities, the independent and dependent variables are nondimensionalised as follows (Chambers et al. 2022 These scalings are chosen to give ∞ 1m ã,t dã = ∞ 1p ã,t dã = 1, such that m ã,t andp ã,t may be considered as probability density functions for the live and apoptotic macrophage populations, respectively. We further define the following dimensionless parameters: The dimensionless model equations are presented in full in Appendix A. Appendix B presents three additional ODEs that describe the time evolution of the total amount of lipid in the system L (t) =Ã M (t) +Ã P (t) +Ñ (t), the average lipid content per live cellĀ M (t) =Ã M (t)

M(t)
, and the average lipid content per apoptotic cellĀ P (t) =Ã P (t) . These equations prove useful for interpreting the simulation outcomes. The numerical method used to solve the model equations is summarised in Appendix C.

Base Case Model Parameterisation
The aim of this study is to investigate how lipid-dependent macrophage behaviour can alter the fate and dynamics of a plaque compared to cases where macrophage behaviour is lipid-independent. Thus, to maintain our focus on the role of lipid-dependence, we vary only the parameter values that characterise the functions g (a) and we fix all parameter values that are not related to these functions (see Table 1). Details on the parameterisation of the functions g (a) are provided at the beginning of Sect. 3. Here, we provide justification for the values of the lipid-independent base case parameters. We estimate that β = ν = 0.05 h −1 . The value for the apoptosis rate β is consistent with estimates from Thon et al. (2018) and Ford et al. (2019b), who fitted models to data from in vitro experiments on macrophage lipid loading. The secondary necrosis rate ν is estimated from observations that cell lysis occurs between 12-24 h after apoptosis (Collins et al. 1997).
Using data from an in vivo study of the monocyte/macrophage inflammatory response in murine myocardial infarction, we estimate γ = 0.01 h −1 and α = 10 4 cells h −1 (Leuschner et al. 2012). The value of the emigration rate γ corresponds to the observation that 5-15% of cells exit the site of inflammation rather than die in situ, where we take our estimate at the upper end of the interval. The value for the maximal macrophage recruitment rate α is a conservative estimate based on the observation that the monocyte/macrophage turnover in the inflamed infarct exceeds 10 6 cells in a 24 h period.
For the efferocytosis rate, we choose η = 2.4×10 −6 cell −1 h −1 , which is consistent with estimates made by Marée et al. (2005) in a study that fits ODE models to in vitro data on macrophage engulfment of apoptotic cells. An accurate estimate of the necrotic lipid ingestion rate θ is not currently available. However, it is believed that the rate of macrophage phagocytosis of apoptotic cells is considerably more efficient than that of necrotic material (Kojima et al. 2017). Hence, we assume that η > θ, and set θ = 1.5 × 10 −7 cell −1 h −1 .
Due to a lack of suitable quantitative data, estimates for the values of the remaining parameters ρ, λ and κ are not currently available. For macrophage proliferation, we consider the two cases ρ = 0 h −1 and ρ = 0.025 h −1 . Here, the larger value ensures compliance with the condition (12), and, in the absence of other factors, predicts a biologically reasonable macrophage population doubling time of approximately 28 h. For λ and κ, we choose dimensionless values that lead to a mild inflammatory response and relatively benign plaque formation in the absence of lipid-dependent effects. By making this choice, we leave scope to investigate whether lipid-dependent macrophage behaviour can, as one might anticipate, promote a more robust inflammatory response and lead to potentially dangerous plaque formation.

Results
The model outlined in Sect. 2 includes lipid-dependent terms for macrophage apoptosis, emigration and proliferation. We stress, however, that it is not our intention to apply the model in its full generality. Rather, we shall consider lipid-dependence in each behaviour individually to investigate how each lipid-dependent behaviour can influence plaque progression. For the functions g (a), we consider a range of parameterisations, both in unscaled and scaled formats. We use a range of parameterisations for the lipid-dependent cell behaviours because there is a lack of appropriate experimental data with which to accurately estimate the relevant quantities. For unscaled simulations, we apply the functions g s (a) or g r (a) exactly as defined in Sect. 2.4. For scaled simulations, we pre-multiply the relevant function by a scaling value that we have found (by simulation) to give G (t) ≈ 1 at steady-state (specifically, we accept G (∞) = 1 ± 0.01). Scaling the functions in this way allows for a consistent comparison of steady-state results from both lipid-dependent and lipid-independent cases because the net (population level) rate of the behaviour of interest is approximately conserved. We focus mainly on steady state results because the time to reach steady state (typically around 100 macrophage lifetimes) is considerably shorter than the lifespan of a plaque. Table 2 (monotonic functions) and Table 3 (non-monotonic functions) summarise the cases that we consider in the following sections. Each table reports function parameterisations, G (∞) values from unscaled simulations and corresponding scaling values required to give G (∞) ≈ 1.

Lipid-Independent Base Case Simulations
Intra-plaque macrophage proliferation has only recently been established as an important contributor to plaque progression (Robbins et al. 2013). The extent of this proliferation is not well characterised but it is understood to vary over the lifetime Table 2 Parameterisations, unscaled G (∞) values and corresponding scaling values for simulations with the monotonic rate modulating function (32) Cell behaviour g (a) Parameter values Asterisks denote that δ γ is divided by the scaling value when the scaling is applied. This ensures that both the scaled and unscaled functions tend towards the same value as a → ∞ Table 3 Parameterisations, unscaled G (∞) values and corresponding scaling values for simulations with the non-monotonic rate modulating function (33) Note that these simulations do not consider macrophage proliferation (i.e. ρ = 0). Asterisks denote that scaling values multiply only the second term on the right-hand side of (33). This ensures that both the scaled and unscaled functions tend towards γ as a → ∞ of a plaque (Lhoták et al. 2016). Given this uncertainty, we shall perform numerical investigations of lipid-dependent cell behaviour both in the absence (Sect. 3.2) and the presence (Sect. 3.3) of macrophage proliferation. To provide a reference point for our lipid-dependent simulations, we first generate base case results where macrophage behaviour is independent of internalised lipid (all g (a) = 1). Steady-state results for these cases, which use only the parameter values in Table 1, are shown in Fig. 1 . For an in-depth understanding of these results, interested readers are referred to the simulations and analysis in Ford et al. (2019a) and Chambers et al. (2022). Here, we provide only a brief summary of some key features of the results.
In the absence of proliferation (ρ = 0), the model recruits a relatively small macrophage population and forms a moderately sized necrotic core (M ≈ 0.24, N ≈ 4.03; Fig. 1a black bars). This reflects our conservative model parameterisation, since we wish to allow for elevated macrophage recruitment in lipid-dependent cases where poorer outcomes are anticipated (e.g. when apoptosis increases with a, or when emigration decreases with a). In the case with proliferation (ρ = 0.5), we observe an almost 3-fold increase in the live cell population (M ≈ 0.70; Fig. 1a (32) with parameter values n β = 2 and a β = 15, δ β = 2 (dashed lines), a β = 12, δ β = 3 (solid lines) or a β = 9, δ β = 4 (dot-dashed lines). Scaling values for the plots in (b) are 0.86, 0.72 and 0.565, respectively bar). This is partly due to proliferation itself, but also due to enhanced recruitment courtesy of the substantial increase in A M (4.70 vs. 1.96). The increased cell population leads to a reduction in both the necrotic core size (N ≈ 1.45) and the average lipid per cell (Ā M =Ā P ≈ 6.7, down from 8.1). The reduction inĀ M due to proliferation can be seen qualitatively by comparing the steady state m (a, t) distributions in Fig. 1b. The case with proliferation (red line) has a considerably larger proportion of cells with very small lipid loads.
For the lipid-independent parameterisation chosen here, there is no single macrophage behaviour that dominates plaque formation. Given that these lipidindependent parameters are fixed for the entirety of the study, we argue that avoiding a single dominant cell behaviour is the most appropriate way to gain a general appreciation of the impact of lipid-dependence on plaque fate and dynamics.

Lipid-Dependent Simulations Without Proliferation
In this section, we neglect macrophage proliferation and investigate, in turn, the role of lipid-dependence in macrophage apoptosis and emigration.

Apoptosis Only
To investigate lipid-dependent apoptosis, we set g γ (a) = 1 and g β (a) = g s β (a). We consider three different unscaled forms for g s β , each of which has δ β > 1 and n β = 2 (Fig. 2 a). Each function reflects an assumption that the likelihood of macrophage apoptosis increases with increasing ingested lipid content (Tabas 2002;Feng et al. 2003). By varying the values of both a β and δ β , we investigate the impact of a mild (a β = 15, δ β = 2), moderate (a β = 12, δ β = 3) or severe (a β = 9, δ β = 4) increase in the apoptosis rate above the reference value with increasing a. (By changing the  t) and N (t) for the lipidindependent base case (g β = 1; black lines) and for three cases with unscaled lipid-dependent apoptosis (g β = g s β ). Lipid-dependent cases have parameter values n β = 2 and a β = 15, δ β = 2 (red lines), a β = 12, δ β = 3 (blue lines) or a β = 9, δ β = 4 (green lines) (Color figure online) parameter values in this way, we anticipate that all scenarios with a β ∈ [9, 15] and δ β ∈ [2, 4] will have solutions that lie within the bounds of those reported below.) Time dependent solutions of the ODE variables for these cases are compared with those for the base case simulation (g β (a) = 1) in Fig. 3. An immediate observation is the interesting behaviour of the case with a β = 9 and δ β = 4. While the other cases all display similar dynamics, this case elicits oscillations that eventually decay towards steady state. Noting that this is the case with the highest net apoptosis rate (dimensionless value G β (t) = 2.367 at steady state; see Table 2), a likely explanation for the oscillations is as follows. Initially, live macrophage numbers M (t) drop to a a 45 and 0 t 1200 in a simulation with unscaled lipid-dependent apoptosis (g β = g s β with n β = 2, a β = 9 and δ β = 4). Note the onset of temporal oscillations in m (a, t) near a = 1 at t ≈ 400. These oscillations reflect repeated bursts of macrophage recruitment that gradually diminish in intensity very low level because cells that ingest even moderate quantities of lipid quickly die to become apoptotic macrophages P (t). The associated conversion of live cell lipid A M (t) to apoptotic cell lipid A P (t) reduces the rate of macrophage recruitment. Low cell numbers and sustained apoptotic lipid generation lead to rapid accumulation of necrotic lipid N (t). Eventually, however, the necrotic lipid pool becomes so vast that the live macrophages begin to ingest lipid at a rate that outstrips the death rate. This allows A M (t) to rise, which stimulates macrophage recruitment and enhances lipid ingestion to shrink the necrotic core. Of course, as more live cells now attain higher quantities of ingested lipid, their apoptosis rates increase further and M (t) drops once again. A new period of growth in N (t) is then initiated until a further (smaller) wave of macrophage recruitment is triggered. This cycle repeats until the magnitude of the oscillations tend to zero. Figure 4presents a corresponding surface plot of the live macrophage distribution m (a, t) in this case. The solution shows temporal oscillations near a = 1, which are associated with the repeated waves of macrophage recruitment.
Long-time solutions of the ODE variables for each simulation indicate that the magnitude and steepness of the increase in g s β correlates with a decrease in M and an increase in each of the other variables (ranging from a marginal rise in A M to a substantial rise in N ). The magnitude and steepness of the change in g s β also appears to be correlated to the time required to reach steady state, which increases from t ≈ 100 in the base case to t ≈ 180 (a β = 15, δ β = 2), t ≈ 350 (a β = 12, δ β = 3) or t ≈ 2000 (a β = 9, δ β = 4). The steady state distributions of live macrophages m (a, t) and apoptotic macrophages p (a, t) for each simulation are presented in Fig. 5 a (note, from Eq. (24), that m (a, t) and p (a, t) differ at steady state only in the case of lipid-dependent apoptosis). For both m and p, we see that increasing the magnitude and steepness of the change in g s β skews the distributions towards larger lipid loads. This seems counter-intuitive, but it appears that the elevated lipid ingestion rate (due to the increase in N and A P ) produces more cells with large lipid loads than are Steady state m (a, t) distributions (left panels) and p (a, t) distributions (right panels) for the lipidindependent base case (g β = 1; black lines) and for three cases with a unscaled or b scaled lipid-dependent apoptosis (g β = g s β ). Lipid-dependent cases have parameter values n β = 2 and a β = 15, δ β = 2 (red lines), a β = 12, δ β = 3 (blue lines) or a β = 9, δ β = 4 (green lines). Primary plots show the results on the interval a ∈ [1, 10] and inset log-log plots show the results on the entire a domain (Color figure online) lost through apoptosis. At steady state, the m and p distributions are related by the expression p (a, ∞) = [see Eq. (24)]. In the simulation with a β = 9 and δ β = 4, this gives a non-monotonic profile for p (a, ∞) with a shallow peak near a = 8.
When interpreting these results, it is important to remember that the net apoptosis rates G β (t) vary considerably between the different simulations (from 1 in the base case to 2.367 at steady state in the case with a β = 9 and δ β = 4). To correct for this difference, we repeat the lipid-dependent simulations with appropriately scaled functions g s β (a) that give G β (t) ≈ 1 at steady state (Fig. 2b). For these simulations, we find that the steady state values of M, P and A M effectively remain fixed, and only A P and N vary across the different cases. (The oscillatory dynamics observed above do not occur in this case, and, while the trend of increasing time to steady state remains, it is much less pronounced than before). Figure 6compares the steady state solutions for A P and N from the base case with those from the scaled lipid-dependent simulations. While the absolute variation in these quantities is much smaller than in Fig. 6 Steady state solutions of the ODE variables A P (t) and N (t) for the lipid-independent base case (g β = 1; black bars) and for three cases with scaled lipid-dependent apoptosis (g β = g s β ). Lipid-dependent cases have parameter values n β = 2 and a β = 15, δ β = 2 (red bars), a β = 12, δ β = 3 (blue bars) or a β = 9, δ β = 4 (green bars). Solutions for M (t), P (t) and A M (t) are omitted from the plot as their values are unchanged across cases (Color figure online) the unscaled scenario, the trend of increasing A P and N with increasing magnitude and steepness of the change in g s β remains. Plots of the corresponding steady state PDE solutions are shown in Fig. 5b. These demonstrate that the trends in the m (a, ∞) and p (a, ∞) distributions with scaled g s β (a) are largely conserved from the unscaled cases (Fig. 5a), although the differences relative to the base case results are less pronounced. Each of the p (a, ∞) distributions presented in Fig. 5b satisfies the relationship p (a, ∞) (a, ∞). Accordingly, the p (a, ∞) distribution in each case corresponds exactly to the distribution of apoptosis events at steady state. Moreover, Eq. (36) shows that, in each case, the steady state average lipid per apoptotic cellĀ P = A P P is given by the steady state G βa value. Given that the steady state P values are essentially identical across all cases, this indicates that the steady state A P values are exactly proportional to the steady state G βa values. The trend observed in Fig. 6 for increasing steady state A P (and N ) with increasing magnitude and steepness of the change in g s β is therefore related to a trend for increasing steady state G βa . The precise reason for this trend of increasing steady state G βa is, however, not entirely clear. It may be correlated to the increase in the (maximum) steepness of g s β , or it may be correlated to the increase in the limiting value of g s β (c.f. Fig. 2b). The fact that the steady state M and A M values are essentially identical across these scaled simulations has interesting consequences when viewed through the lens of Eqs. (34) and (35). For the scenarios simulated here, Eq. (34) reduces to: The total system lipid L (t) has explicit dependence only on M (t) and A M (t). Interestingly, however, our simulation results indicate that the steady state L value is not uniquely determined by the steady state M and A M values (i.e. even though steady state M and A M are essentially fixed across simulations, there are differences in steady state L due to the variation in steady state A P and N ). This observation demonstrates that, in terms of total system lipid, the model plaque carries the weight of its history. for the lipid-independent base case (g β = 1; black lines) and for three cases with scaled lipid-dependent apoptosis (g β = g s β ). Lipid-dependent cases have parameter values n β = 2 and a β = 15, δ β = 2 (red lines), a β = 12, δ β = 3 (blue lines) or a β = 9, δ β = 4 (green lines). Although all simulations reach similar steady state M and A M values, the paths taken to get there vary and this has consequences for the total quantity of lipid retained in the system (Color figure online) That is, the observed increase in steady state L with decreasing a β and increasing δ β can be explained by the historical accumulation of lipid due to periods of reduced lipid removal by emigration (i.e. reduced A M (t)) or increased lipid addition by recruitment (i.e. increased F (t) or, equivalently, increased A M (t) − M (t)). The time courses of the relevant M (t) and A M (t) solutions (Fig. 7) suggest that reduced lipid removal by emigration is the predominant mechanism in this case.
Unlike L (t), the steady state solution forĀ M (t) remains fixed across all four simulations. In each of these cases, the steady state solution of Eq. (35) can be expressed as: where all time-dependent variables take their steady state values. As steady state M andĀ M have fixed values for all simulations, only terms two, three and six on the left-hand side of (22) change in each case. Thus, steady stateĀ M is unchanged across all scaled simulations because any reduction in average lipid per cell due to lipiddependent apoptosis (term six) is exactly compensated by an increase in necrotic and efferocytic lipid consumption (terms two and three, respectively). This finding supports our earlier interpretation of the unscaled simulation results, where we found that the direct contribution of lipid-dependent apoptosis to reducingĀ M (t) was counteracted by other factors.

Emigration Only
To study lipid-dependent emigration, we set g β (a) = 1 and consider that g γ (a) takes either the monotonic form g s γ (a) with δ γ < 1, or the non-monotonic form g r γ (a). In the monotonic case, we assume that the rate at which macrophages leave the plaque decreases with increasing lipid load (  In the non-monotonic case, we retain this argument. However, we additionally assume a substantially reduced emigration rate for macrophages with little or no accumulated lipid. This assumption reflects a range of considerations, including the innate propensity for macrophages to pursue foreign bodies, and the fact that emigrating macrophages typically traverse the plaque before exiting to the media (Llodrá et al. 2004). Either way, we anticipate that newly-recruited macrophages are unlikely to leave the plaque without first ingesting at least some lipid. Combined with the assumption of reduced emigration for macrophages with large lipid loads, this leads us to consider that the lipid-dependent emigration rate may peak at some intermediate lipid quantity. Note, however, that the argument regarding low rates of emigration for macrophages with small lipid loads may not hold if the lipid load can be reduced by proliferation. As such, we will only consider the function g r γ (a) for the model without macrophage proliferation.
We first investigate monotonic lipid-dependence in macrophage emigration and consider three alternative unscaled forms for g s γ (a) (Fig. 8a). Each function has δ γ = 0.1 and n γ = 1.5, while we vary a γ to study the impact of a gentle (a γ = 24), moderate (a γ = 18) or steep (a γ = 12) reduction in the emigration rate as a increases over low values. We deliberately choose a small n γ value and non-zero δ γ to maintain feasibility in the numerical solution of the equations. When the functions g s γ decline rapidly and/or tend towards zero as a → ∞, a non-negligible proportion of cells can acquire extremely large lipid loads (a 1000). For numerical accuracy, such cases need to be solved on extended domains at significant computational expense. We avoid such scenarios here, as this increase in computational effort is unlikely to produce additional physical insight. While physically we may expect effectively zero emigration of macrophages with very large lipid loads, our chosen parameterisation ensures that cells are still fifty times more likely to die than to emigrate as a → ∞. (c) (d)

Fig. 9 Steady state solutions of the ODE variables M (t), P (t), N (t), A P (t) and A M (t)
for the lipidindependent base case (g γ = 1; black bars) and for several scenarios with lipid-dependent emigration. Lipid-dependent simulations use a unscaled or b scaled monotonic lipid-dependence (g γ = g s γ ) and c unscaled or d scaled non-monotonic lipid-dependence (g γ = g r γ ). Monotonic lipid-dependent cases have parameter values n γ = 1.5, δ γ = 0.1 and a γ = 12 (red bars), a γ = 18 (blue bars) or a γ = 24 (green bars). Non-monotonic lipid-dependent cases have parameter values γ = 0.1, k γ = 1, q γ = 2 and b γ = 3 (red bars), b γ = 6 (blue bars) or b γ = 9 (green bars) (Color figure online) Steady state results for the ODE variables in each lipid-dependent case are compared with those from the base case (g γ = 1) in Fig. 9 a (note the use of several vertical axis scales for ease of visualisation). These results show all five variables trending upwards with decreasing a γ (or, equivalently, with decreasing net emigration rate G γ (∞)). These upward trends are of varying degrees, with the case for a γ = 12 showing approximately 1.3-, 1.8-, 3.1-, 5.5-and 13.7-fold increases in P, N , M, A P and A M , respectively, versus the base case. The increase in steady state M in each lipid-dependent case is partly due to reduced net emigration, but mainly due to increased cell recruitment as a consequence of the substantial rise in A M . The (c) (d) Fig. 10 Steady state m (a, t) distributions for the lipid-independent base case (g γ = 1; black lines) and for several scenarios with lipid-dependent emigration. Lipid-dependent simulations use a unscaled or b scaled monotonic lipid-dependence (g γ = g s γ ) and c unscaled or d scaled non-monotonic lipid-dependence (g γ = g r γ ). Monotonic lipid-dependent cases have parameter values n γ = 1.5, δ γ = 0.1 and a γ = 12 (red lines), a γ = 18 (blue lines) or a γ = 24 (green lines). Non-monotonic lipid-dependent cases have parameter values γ = 0.1, k γ = 1, q γ = 2 and b γ = 3 (red lines), b γ = 6 (blue lines) or b γ = 9 (green lines). Primary plots show the results on the interval a ∈ [1, 10] and inset log-log plots show the results on the entire a domain (Color figure online) steady state A M values increase so significantly because, as cells accumulate more and more lipid, they become less likely to emigrate. This is shown in the inset of Fig. 10a by the (relatively) large proportions of cells with very large lipid loads in the steady state m (a, t) distributions. As these heavily lipid-loaded cells are much more likely to undergo apoptosis than to emigrate, A P is increased and this ultimately fuels an increase in N . These enlarged necrotic lipid pools emerge despite the significant increase in overall lipid consumption afforded by the increase in M.
To correct for the variation in the net steady state emigration rates in the above lipid-dependent cases, we repeat our study with appropriately scaled g s γ (a) functions. The functions are scaled in such a way that the limiting value of each g s γ remains (a) (b) Fig. 11 a Unscaled and b scaled rate modulating functions for macrophage emigration g γ (a) = g r γ (a). Plots correspond to Eq. (33) with parameter values γ = 0.1, k γ = 1, q γ = 2 and b γ = 3 (dashed lines), b γ = 6 (solid lines) or b γ = 9 (dot-dashed lines). Scaling values for the plots in (b) are 1.9, 1.99 and 2.08, respectively. These scaling values multiply only the second term in Eq. (33) so that each g r γ retains the limiting value of the unscaled functions at 0.1. In practice, this involves dividing δ γ by the scaling value when the scaling is applied. We make this assumption because the model results are particularly sensitive to the limiting value of g s γ and we wish to retain consistency for comparison with the unscaled cases. Steady state solutions of the ODE variables are presented in Fig. 9b. Unlike lipid-dependent apoptosis, where the steady state values of M, P and A M were all unchanged from their base case values upon scaling, here we see that all five variables are different from the base case. Indeed, the trends in the ODE solutions are identical to those in the unscaled cases, although with slightly smaller variation from the base case values. This observation correlates with the steady state m (a, t) distributions in Fig. 10b, which become slightly closer to the base case distribution than those in the unscaled case (Fig. 10a). These results with scaled g s γ demonstrate that the particular form of the lipid-dependent emigration function (parameterised here by a γ and an appropriate scaling value) can have considerable influence on the long-term plaque composition.
We now consider the case of non-monotonic lipid-dependent emigration, where the function g r γ (a) encodes a reduced emigration rate for macrophages with small lipid loads. We fix the parameter values k γ = 1, q γ = 2 and γ = 0.1 (note the consistency in the limiting value of g γ ), and vary b γ such that g r γ has its peak at a = 4 (b γ = 3), a = 7 (b γ = 6) or a = 10 (b γ = 9). Plots of these three unscaled g r γ are shown in Fig. 11 a, where we see that each rightward shift in the peak reduces the rate of decline of the function with increasing a. The steady state ODE solutions and m (a, t) distributions generated with these functions are shown in Figs. 9c and 10c, respectively. An immediate observation is that the plot of the ODE results is visually very similar to that for the unscaled monotonic case (Fig. 9a). This, however, is coincidental and should not be regarded as significant. Indeed, the lipid-dependent case with the best outcome (smallest N ) in Fig.9a is the one with the highest net emigration rate at steady state (a γ = 24, G γ (∞) ≈ 0.83), whereas the lipid-dependent case with the Fig. 12 Plot comparing how emigration events g γ (a) m (a, ∞) (solid lines) and associated quantities of removed lipid g γ (a) am (a, ∞) (dashed lines) are distributed with respect to a at steady state in simulations with monotonic lipid-dependent emigration g γ = g s γ (a γ = 24; red lines) and non-monotonic lipiddependent emigration g γ = g r γ (b γ = 9; blue lines). The net emigration rate for the non-monotonic case is considerably smaller than that for the monotonic case (G γ = ∞ 1 g γ m da ≈ 0.50 vs. 0.83), but the increased proportion of emigration events with larger lipid loads produces a similar net lipid removal rate (G γ a = ∞ 1 g γ am da ≈ 6.29 vs. 6.61). The average lipid removed per cell is approximately 8.0 for the monotonic case and approximately 12.6 for the non-monotonic case best outcome in Fig. 9c is the one with the lowest net emigration rate at steady state (b γ = 9, G γ (∞) ≈ 0.50). This observation demonstrates that when emigration events are skewed towards macrophages with larger lipid loads, the efficiency of lipid removal from the system can be substantially improved (Fig. 12). The case with b γ = 9 has a smaller steady state N value than the base case (where G γ = 1), which presumably reflects a preferential balance in the steady state M and A P values. However, the exact mechanism that underlies this result is difficult to ascertain.
Given that non-monotonic lipid-dependent emigration can remove substantial amounts of lipid from the plaque even at reduced net emigration rates, it seems reasonable to expect an improvement in outcomes when the net emigration rate is scaled to the reference value. (See Fig. 11b for plots of the scaled lipid-dependent functions g r γ . We again preserve the limiting value of 0.1 by multiplying only the non-constant part of (33) by the scaling value.) This is indeed true in some cases, but the steady state ODE results in Fig. 9d portray a greater subtlety in the model's response to this scaling. Relative to their corresponding unscaled cases (Fig. 9c), each scaled case shows reduced steady state M, P, A M and A P values. Steady state N also decreases for b γ = 3 and b γ = 6 (here falling below the base case value), but for b γ = 9 the core size increases. It seems that, in this case, the extent of lipid removal from the plaque is so substantial that recruitment is reduced and too few macrophages remain in the plaque to resolve the necrotic core. This explanation is corroborated by the steady state m (a, t) distribution in Fig. 10d (green line). Compared to the corresponding result for the unscaled case (Fig. 10c), we see large reductions in both the proportion of cells with very large lipid loads (see inset) and the proportion of cells with very small lipid loads (a ≈ 1). The cases for b γ = 3 and b γ = 6 (red and blue lines, respectively) show similarly large reductions in m at the top end of their distributions, but at the lower end both distributions increase.

Lipid-Dependent Simulations with Proliferation
In this section, we introduce macrophage proliferation into the system by setting ρ = 0.5. First we assume that proliferation is independent of lipid load. We briefly revisit the scenarios from Sect. 3.2 to study the manner in which proliferation (with its tendency to reduce average lipid loads) interacts with lipid-dependence in the other cell behaviours. We then revert to lipid-independent apoptosis and emigration, and conclude the results by investigating the impact of lipid-dependence in proliferation itself.

Apoptosis Only or Emigration Only
We set g ρ (a) = 1, and select one lipid-dependent case from each of Sects. 3.2.1 and 3.2.2. These are g β = g s β with a β = 12, δ β = 3 and g γ = g s γ with a γ = 18, respectively. As before, we allow lipid-dependence in only one cell behaviour at a time, and we investigate each lipid-dependent case using both scaled and unscaled functions. (Note that g β and g γ generally require new scaling values to account for the impact of proliferation; see Table 2). Here, we shall primarily focus on the scaled cases since preserving the net steady state apoptosis or emigration rate facilitates comparison against the earlier scaled cases without proliferation. The unscaled cases do, however, provide some interest and we comment on these below.
For the unscaled lipid-dependent apoptosis case, Table 2 shows that the inclusion of proliferation acts to reduce the net steady state apoptosis rate G β by approximately 8% from 1.576 to 1.445. This reflects the fact that proliferation acts to reduce individual cell lipid loads, thereby reducing the likelihood of apoptosis (i.e. daughter cell apoptosis rates g s β (a) are always less than parent cell apoptosis rates g s β (2a − 1)). Note, however, that a tangible reduction in the daughter cell apoptosis rate requires g s β (a) g s β (2a − 1), and this occurs only for daughter cells with a sufficiently small. The unscaled lipid-dependent emigration case produces the rather curious result that proliferation actually reduces the net steady state emigration rate (G γ (∞) drops from 0.7812 to 0.7650; see Table 2). This is counter-intuitive because the associated reduction of lipid loads would be expected to increase the likelihood of emigration (i.e. g s γ (a) > g s γ (2a − 1) for all a). This quirk appears to arise due to the proliferation of cells with very large lipid loads. Proliferation of such cells fails to significantly elevate daughter cell emigration rates because g s γ (a) remains similar to g s γ (2a − 1) for large a. Thus, lipid loaded cells with very low emigration rates perpetuate in the system and skew the G γ value upwards. Despite this unexpected increase in G γ , proliferation does, as expected, reduce the average lipid per cell at steady state.
Returning to the scaled lipid-dependent functions, Fig. 13 presents the steady state ODE solutions for the lipid-dependent and base case simulations both in the absence (a) (b)

Fig. 13 Steady state solutions of the ODE variables M (t), P (t), N (t), A P (t) and
A M (t) for lipidindependent (g β = g γ = 1; black bars) and scaled lipid-dependent cases a without proliferation (ρ = 0) and b with proliferation (ρ = 0.5). Red bars show results with lipid-dependent apoptosis (g β = g s β , g γ = 1) using parameter values n β = 2, δ β = 3 and a β = 12. Blue bars show results with lipid-dependent emigration (g γ = g s γ , g β = 1) using parameter values n γ = 1.5, δ γ = 0.1 and a γ = 18 (Color figure online) (Fig. 13a) and presence (Fig. 13b) of macrophage proliferation. Although the two solution sets are quantitatively different, the lipid-dependent results display the same qualitative trends regardless of the particular proliferation rate. Lipid-dependent apoptosis increases N and A P (M, P and A M remain fixed), while lipid-dependent emigration increases the solutions of all five ODE variables. Figure 14 shows the corresponding steady-state m (a, t) distributions for the cases without proliferation (Fig. 14a) and with proliferation (Fig. 14b). As proliferation adds a second non-local effect into the model, the precise features of the plots in Fig. 14b are difficult to interpret. However, as expected based on Eq. (35), proliferation tends to shift the m (a, t) distributions towards lower accumulated lipid loads. For lipid-dependent apoptosis, proliferation produces a relatively large increase in the proportion of cells with small lipid loads. However, for large a, the impact of proliferation seems relatively minimal. This likely reflects the substantial increase in the apoptosis rate with increasing a, which effectively acts to suppress proliferation events.

Proliferation Only
Finally, we investigate the role of lipid-dependent macrophage proliferation. In this section, we set g ρ (a) = g s ρ (a) and fix g β (a) = g γ (a) = 1. Once again, we consider three different forms for the unscaled lipid dependent functions (Fig. 15 a). These functions have n ρ = 2, δ ρ = 0 and a ρ = 4, 9 or 14. By varying a ρ , we study the effect of different rates of decline in the macrophage proliferation rate with increased lipid loading (Kim et al. 2018). The steady state ODE solutions generated using these functions are presented in Fig. 16. Here, we observe marginal decreases in P and A P , marked decreases in M and A M , and a marked increase in N with decreasing a ρ . These variations in the steady state solutions across the lipid-dependent cases are not explicitly related to the particular form of g ρ , but simply reflect the associated changes in the net proliferation rates ρG ρ (∞). The case with a ρ = 4, for example, has net steady state proliferation rate ρG ρ ≈ (0.5) (0.6377) ≈ 0.319. A lipid-independent simulation (G ρ = 1) with proliferation rate ρ = 0.319 would therefore give almost identical steady state ODE

Fig. 16 Steady state solutions of the ODE variables M (t), P (t), N (t), A P (t) and
A M (t) for the lipidindependent base case (g ρ = 1; black bars) and for three cases with unscaled lipid-dependent proliferation (g ρ = g s ρ ). Lipid-dependent cases have parameter values n ρ = 2, δ ρ = 0 and a ρ = 4 (red bars), a ρ = 9 (blue bars) or a ρ = 14 (green bars) (Color figure online) solutions, albeit with different temporal dynamics and different steady state PDE solutions.
It is interesting to observe that the G ρ values in these lipid-dependent simulations remain fairly close to 1, even when g s ρ declines rapidly with a. This suggests that the model is relatively insensitive to the value of the parameter a ρ , provided it is not too close to 1. This makes intuitive sense because proliferation generally acts to keep cell lipid loads close to a = 1, and, provided a ρ is not too small, the lipid-dependent proliferation rate remains high throughout this region (i.e. g s ρ ≈ 1). Contrastingly, while lipid-independent proliferation may allow some cells to acquire small lipid loads from large lipid loads through multiple proliferation events, this phenomenon will be suppressed in the lipid-dependent cases due to the decline in g s ρ for large a. Consistent with their definition, simulations using the scaled g s ρ (Fig. 15b) lead to net proliferation rates approximately equal to ρ at steady state. Accordingly, each of these simulations produces steady state ODE results that are indistinguishable from the base case. We therefore omit these ODE results and focus instead on the corresponding steady state m (a, t) distributions (Fig. 17 ). We also omit the m (a, t) distributions generated using the unscaled g s ρ . These differ slightly from the scaled g s ρ results for small a, but the distributions become increasingly similar as a increases. Figure 17 shows that there are few discernible differences between the base case steady state m (a, t) distribution and the distributions generated using the scaled g s ρ . One exception is that the lipid-dependent cases all have a larger proportion of cells with lipid load a > 50. This is because these cells have substantially reduced proliferation rates and cannot readily divide to reduce their lipid load. In the case with a ρ = 4, there is a subtle change in the form of the distribution with the emergence of a small peak near Fig. 17 Steady state m (a, t) distributions for the lipid-independent base case (g ρ = 1; black lines) and for three cases with scaled lipid-dependent proliferation (g ρ = g s ρ ). Lipid-dependent cases have parameter values n ρ = 2, δ ρ = 0 and a ρ = 4 (red lines), a ρ = 9 (blue lines) or a ρ = 14 (green lines). The primary plot shows the results on the interval a ∈ [1, 10] and the inset log-log plot shows the results on the entire a domain. All four distributions coincide at a = 1 because the values of the ODE variables in the boundary condition (30) are almost identical in each case (Color figure online) a = 1. This is likely due to the locally elevated proliferation rate, which alters the relative balance between proliferation and efferocytosis (see Fig. 7 in Chambers et al. (2022)). Overall, we conclude from the results in this section that the outcome of plaque formation is probably more sensitive to the net population-level proliferation rate than to the particular lipid-dependent distribution of proliferation rates across the population.

Discussion
Lipid consumption is known to alter macrophage behaviour but the implications for atherosclerotic plaque formation are not well understood. In this paper, we establish a novel modelling framework to study these implications. We develop a structured population model of plaque macrophage lipid accumulation in which macrophages are classified by their internalised lipid load a and behave in a lipid-dependent manner. This work builds upon the partial integro-differential equation models recently developed in Ford et al. (2019a) and Chambers et al. (2022). As in these earlier works, we model how live macrophage behaviour feeds into apoptotic macrophage population dynamics and necrotic core formation through mechanisms such as efferocytosis, post-apoptotic necrosis and necrotic core consumption. A key difference in the current model is that the lipid-averaged ODE subsystem cannot be easily decoupled from the governing PDEs.
We focus our investigation on how lipid-dependent macrophage apoptosis, emigration and proliferation can influence plaque progression. This work is based on results that show dysfunctionality in macrophages that have large ingested lipid loads (Tabas 2002;Moore et al. 2013;Yin and Heit 2021). Consequently, we assume that increasing lipid loads cause increased rates of apoptosis and reduced rates of both emigration and proliferation. We do not consider lipid-induced dysfunction in macrophage phagocytic or efferocytic capacity. Thus, heavily lipid-loaded macrophages in our simulations consume extracellular lipid and engulf apoptotic bodies at the same rate as cells with small ingested lipid loads. This is a limitation of the current study that we will address in future by including lipid-dependent phagocytosis and efferocytosis terms in the model. We anticipate that the inclusion of these terms will have dual benefit. Not only will the model become more realistic, but lipid consumption rates that decrease with increasing lipid load will reduce the domain sizes required for accurate numerical solutions. Model simulations will therefore be less computationally demanding and numerical techniques such as the non-uniform gridding strategy implemented here may no longer be required.
Lipid-dependent macrophage behaviour is encoded in the model through the dimensionless functions g β (a), g γ (a) and g ρ (a), which modulate the reference rates of apoptosis, emigration and proliferation, respectively. These functions typically take a monotonic form, but for emigration we also consider a non-monotonic lipiddependence. In the absence of appropriate data to parameterise these functions, we simulate a range of scenarios with different rates of approach to their (finite) limiting values. This is a reasonable strategy as the lipid-dependence of each behaviour in vivo is likely to be controlled by several different factors including the activation or disruption of unique signalling pathways (Feng et al. 2003;van Gils et al. 2012;Robbins et al. 2013). Thus, the rate of each individual macrophage behaviour may be altered in different ways, at different times and at different internalised lipid loads. In our analysis, we have made the simplifying assumption that only one of apoptosis, emigration or proliferation may be lipid-dependent at any one time. While this may be unrealistic in practice, the benefit of this assumption is that we can study the influence of each lipid-dependent behaviour on plaque progression. Of course, simulations with more than one lipid-dependent behaviour can be easily performed within our model framework. However, the increased difficulty in interpretation of such results limits the insight that can be gained.
For each scenario that we model in this paper, we fix the lipid-independent parameter values and perform several simulations using both unscaled and scaled functions g (a). Unscaled simulations use the exact functional forms (32) or (33) with an appropriate set of lipid-dependent parameter values. Each corresponding scaled simulation uses the same lipid-dependent parameter values, but the function g is multiplicatively scaled such that the net rate of the behaviour of interest matches the reference rate at steady state (mathematically, this is expressed as G (∞) = ∞ 1 g (a) m (a, ∞) da = 1 where g (a) = 1). Each approach has unique advantages and disadvantages. Unscaled simulations are relevant for understanding the impact of variability in lipid-dependent macrophage behaviour and/or the response of the model system to interventions that alter the lipid-dependence. However, as the net steady state rate of the behaviour of interest is not conserved from case to case, it is difficult to interpret the impact of the lipid-dependence itself. By correcting for this, scaled simulations provide a more consistent and insightful means of comparison across simulated scenarios. However, the scaling value required in a given case is not known a priori and must be determined by a process of estimation and refinement through repeated simulation. As such, the scaling process itself can be computationally demanding.
Our results and analysis indicate some interesting population level differences in the relative influence of lipid-dependent apoptosis, emigration and proliferation on plaque progression. For unscaled lipid-dependent proliferation simulations with net steady state proliferation rate ρG ρ (∞), we find that the steady state ODE solutions M, P, A M , A P and N can always be matched by an equivalent lipid-independent simulation with the constant proliferation rate ρG ρ (∞). Consequently, for all scaled lipid-dependent apoptosis simulations (G ρ (∞) ≈ 1), we obtain essentially identical steady state ODE results. This phenomenon is unique to lipid-dependent proliferation and equivalent outcomes are not seen for lipid-dependent apoptosis or lipid-dependent emigration simulations. Mathematically, this reflects the presence in the model equations of the terms G βa and G γ a . These terms appear because the net rate of lipid removal from the live cell population by lipid-dependent apoptosis or lipid-dependent emigration depends on the distribution of lipid across the population. There is no equivalent G ρa term in the model equations because the quantity of lipid added to the system by a proliferation event (a 0 in dimensional terms) is always fixed and does not depend on the internalised lipid content of the parent cell. Thus, compared to lipid-dependent proliferation, the model equations for lipid-dependent apoptosis or lipid-dependent emigration have additional nonlinear terms that increase the complexity of the underlying dynamics. This is emphasised by Eq. (35), where we see that lipid-dependent apoptosis and lipid-dependent emigration can dynamically alter A M (t) in a way that is not observed at all in lipid-independent cases.
In our unscaled lipid-dependent apoptosis simulations, we observe the emergence of decaying oscillations in the ODE variables when the increase in g β (a) is sufficiently large in both magnitude and steepness. Oscillatory solutions of this type appear to be unique to the lipid-dependent model as they have not been observed in earlier lipid-independent approaches (Ford et al. 2019a;Chambers et al. 2022). Two factors that appear to be necessary to initiate these oscillations are: (1) macrophages ingest (primarily necrotic) lipid at such a rate that A M grows irrespective of a concurrent increase in the net apoptosis rate; and (2) this increase in A M stimulates an increase in macrophage recruitment. The latter requirement may only be satisfied for certain values of the dimensionless parameter κ. We use κ = 5 in our simulations. However, for a smaller κ value, the recruitment rate would probably be less sensitive to changes in A M and it is possible that oscillations may not emerge at all.
The observation that oscillatory solutions require a sufficiently rapid rate of lipid ingestion relative to apoptosis raises an interesting question about what might happen if g β does not saturate and instead approaches infinity at some finite a; that is, if there exists a certain lipid load above which macrophages become unviable. It is not clear whether oscillations could occur in this case because A M may be unable to grow due to (potentially) very rapid lipid-induced apoptosis. Non-saturating lipid-dependent death rates may be generally interesting to investigate in this model because, together with lipid ingestion rates that decrease with a, they offer another means to reduce (and, indeed, cap) the length of the a domain.
Oscillatory solutions have also been observed in a plaque formation model by Bulelzai et al. (2014). While it is possible that there may be underlying similarity in the mathematical structure of the models, the biological mechanisms that produce these oscillations in the respective models are quite distinct. These independent observations naturally raise the question as to whether such oscillations are biologically realistic. To our knowledge, oscillatory behaviour in real plaques has never been reported experimentally. However, most plaque data is collected ex vivo and in vivo plaques are typically not observed with sufficient resolution or sufficiently often for such an observation to be made. The oscillations reported here reflect a particular set of modelling assumptions, and, in practice, there may be other factors (e.g., spatial effects, other lipid-dependent behaviours) that act to modify the dynamics. However, the dynamic feedback inherent in monocyte/macrophage recruitment, as well as the multiple timescales inherent in plaque progression, do support the plausibility of oscillatory dynamics in real plaques.
For scaled lipid-dependent apoptosis, we find the surprising result that the steady state M, P and A M values remain unchanged across all simulated cases. Our analysis shows that this happens because, at steady state in each simulation, the net rate of lipid loss from live cells due to apoptosis (G βa ) is exactly balanced by the rate of lipid gain via efferocytosis and necrotic lipid consumption (η A P + θ N ). The balance between these terms fixes the steady state A M value in each simulation, and fixed steady state M and P follow because the recruitment rate, net apoptosis rate and all other relevant rates remain unchanged. Of course, while the relationship between G βa and η A P + θ N causes steady state M, P and A M to remain fixed, the very same relationship causes steady state A P and N to vary. The steady state values of G βa , A P and N all increase with increasing magnitude and steepness of change in scaled g s β . These results suggest that plaques whose cells have increased susceptibility to the cytotoxic effects of lipid accumulation can have larger necrotic cores at steady state even if the net macrophage apoptosis rate is held fixed. Based on the ODE for total system lipid L (t), we further propose that this result reflects excess lipid accumulation in the plaque over time, primarily due to a sustained period of reduced lipid removal by emigration.
Our results for lipid-dependent emigration are arguably the most difficult to analyse and interpret because, even in the case of scaled g γ (a), none of the steady state ODE solutions remain fixed across our simulated cases. This reflects the influence of the term G γ a on the dynamics of A M , which then affects the other ODE variables via the macrophage recruitment rate. Although not considered here, our results suggest that it may be worthwhile to perform a complementary set of lipid-dependent simulations where the functions g γ are scaled such that, at steady state, the quantity G γ a matches its base case value; that is, scale the g γ to match the rate of lipid removal by macrophage emigration at steady state rather than the rate of macrophage emigration itself. While this approach is unlikely to keep any of the individual steady state ODE solutions at their base case values, it may give additional insight into the intricacies of the simulated outcomes.
In our simulations with unscaled monotonic lipid-dependent emigration, we find that the steady state values of the ODE variables all increase with increasing rate of decline of g γ (decreasing a γ ). This is unsurprising because, as a γ decreases, the net emigration rate also decreases and the system therefore retains more cells and more lipids. In simulations with scaled g γ , we observe the exact same trend in the steady state ODE values. This shows that the reduced net emigration rate in the unscaled cases is only partially responsible for the observed outcomes and, in fact, it is the form of the lipid-dependence that contributes much of the increase in the steady state ODE values relative to the base case. This observation reflects the fact that any cell that fails to emigrate for a small becomes increasingly unlikely to emigrate at all. Cells can therefore persist in the system, acquiring lipid and eventually undergoing apoptosis. This leads to large ingested lipid quantities in both the live and apoptotic cell populations, which ultimately feed into enlarged necrotic cores. Our results show that the detrimental effects of monotonically decreasing lipid-dependent emigration can be substantial. Scaled simulations predict up to a 40% increase in necrotic core size despite a 2-to 3-fold increase in the net rate of necrotic core consumption.
As well as monotonic lipid-dependent emigration, we consider non-monotonic g γ where macrophages with moderate lipid loads have the highest emigration rates. These g γ impose reduced rates of emigration for macrophages with very small lipid loads on the basis that such cells (in the absence of proliferation) seem unlikely to exit the plaque (Llodrá et al. 2004). Our results show that non-monotonic lipid-dependent emigration can be beneficial relative to the monotonic lipid-dependence. As cells are more likely to emigrate with a relatively large lipid load, this not only removes more lipid from the system but it also (in scaled cases at least) reduces the likelihood that cells acquire large lipid loads before dying. Our simulations with parameter value b γ = 9 are particularly interesting. In the unscaled case, we see a reduction in steady state N relative to the base case. This is particularly remarkable because: (1) the steady state G γ is only half of the base case value; and (2) the emigration rate of cells with very large a is up to 10 times smaller than in the base case. In the equivalent scaled case, however, we see a larger steady state N relative to the base case. This is because extensive lipid removal from the system by emigration inhibits the immune response. Taken together, these results suggest that the presence of at least some heavily lipid-loaded (and highly inflammatory) macrophages in the plaque may be beneficial for stimulating the recruitment that is required for necrotic core consumption. Observations such as this may be useful in the interpretation of results from experimental studies on plaque resolution and regression (Rahman and Fisher 2018).
Our lipid-dependent emigration results provide new insight into how the distribution of emigration events can influence the extent of lipid removal from the plaque and the overall plaque progression. These results, however, should be interpreted with caution. We assume that the likelihood of macrophage emigration depends only on lipid load, whereas, in practice, this relationship would be time-and space-dependent due to factors such as plaque size and the strength of emigratory signals relative to other migratory cues (e.g. "find me" signals from apoptotic cells (Kojima et al. 2017)). The inclusion of spatial effects in future models will allow a better understanding of how macrophage emigration events may be distributed with respect to internalised lipid loads, both over time and at steady state.
In the current work, we perform two sets of simulations that include macrophage proliferation. We first investigate how lipid-independent proliferation interacts with lipid-dependence in the other kinetic terms, and we then investigate the impact of a monotonically decreasing lipid-dependent proliferation rate. When lipid-independent proliferation is combined with scaled lipid-dependent apoptosis or emigration, we observe, as expected, a general reduction in necrotic core size and average cell lipid loads. Overall, however, the qualitative trends in the steady state ODE solutions remain the same as in the proliferation-free cases. When g β is unscaled, we make the interesting observation that proliferation can improve the outcome of plaque formation by reducing the net apoptosis rate. Although the observed reduction at steady state is relatively small (around 8%), the lipid-dependent apoptosis rates range from 2-6 times larger than the proliferation rate (and, as observed in proliferation-free cases, lipid-dependent apoptosis tends to drive cells to higher a by increasing lipid consumption). For simulations with lipid-dependent proliferation, we identify two main findings. First, that the steady state ODE solutions depend only on the net steady state proliferation rate ρG ρ , and, second, that the steady state distribution of lipid across the live cell population is relatively insensitive to the form of g ρ (provided steady state G ρ is not too far from 1).
Although we consider a monotonically decreasing form for g ρ , there exists experimental evidence that the rate of lipid-dependent macrophage proliferation may peak at an intermediate lipid content (Xu et al. 2015). This form of lipid-dependent proliferation can be easily included in the model by allowing g ρ to take the non-monotonic form defined in Eq. (33). An interesting implication of this non-monotonic lipid-dependence is that, unlike monotonically decreasing lipid-dependent proliferation, it can naturally explain the observation that the contribution of plaque macrophage proliferation increases over time (Lhoták et al. (2016); i.e. as plaque-resident macrophages gradually accumulate lipid, the net population proliferation rate would grow). Of course, many other factors may contribute to this phenomenon, including the presence of a fibrous cap in mature plaques that may inhibit the rate of monocyte recruitment from the bloodstream. We have performed preliminary model investigations using nonmonotonic lipid-dependent proliferation, but we omit the results here because we find that they are not substantially different to those with monotonic lipid-dependence. However, in simulations that consider more than one lipid-dependent behaviour at a time (e.g. lipid-dependent proliferation and lipid-dependent apoptosis), we envisage that the particular form of lipid-dependent proliferation may have a more significant impact on the outcome of plaque progression. For example, a lipid-dependent macrophage proliferation rate that peaks close to the lipid load at which the lipiddependent apoptosis rate begins to rise would potentially optimise the protective effect of proliferation against the otherwise detrimental effects of lipid cytotoxicity.
Before concluding, we emphasise that the aim of this study is to use mathematical modelling to identify ways in which lipid-dependent macrophage behaviour may modify atherosclerotic plaque fate and dynamics. We also wish to understand, from a mathematical perspective, the mechanisms that underlie these modifications. By design, we consider the implications of several different lipid-dependent macrophage behaviours in isolation. This approach provides an ideal way to achieve the stated aims. However, as this model setup is somewhat artificial, we remain naturally cautious when interpreting our findings with respect to in vivo plaque formation. Real plaques are likely to have multiple different lipid-dependent macrophage behaviours at play at any given time and, hence, direct extrapolation of our findings to existing in vivo observations is challenging. The real significance of this work is that it demonstrates that the effects of lipid-dependent macrophage behaviour on plaque formation can be substantial, and therefore highlights that more detailed experimental measurements (including, for example, macrophage internalised lipid distributions) will almost certainly be required to fully comprehend the mechanisms of in vivo plaque progression.

Conclusions
This paper extends the work of Ford et al. (2019a) and Chambers et al. (2022) by developing a lipid-structured model of atherosclerotic plaque macrophages in which the rates of macrophage apoptosis, emigration and proliferation are modulated by the internalised lipid load. We model these lipid-dependent behaviours using dimensionless modulating functions whose features align qualitatively with a range of experimental observations. Our results indicate, particularly for apoptosis and emigration, that variations in macrophage behaviour across lipid loads can substantially alter plaque fate relative to cases without lipid-dependence. We find that these changes to plaque fate are difficult to predict because the lipid-dependent behaviours introduce subtle, nonlinear effects that feed back into the model through phagocytosis, efferocytosis and recruitment. This work provides new insight into how macrophage lipid accumulation can shape atherosclerotic plaque progression and highlights the importance of mathematical modelling as a tool in the scientific study of atherosclerotic plaque macrophages. and the initial conditions are: We set a σ = 0.5 for all simulations in this paper. The functions that control lipid-dependent cell behaviour are: g s (a) = (a − 1) n + δ (a − 1) n (a − 1) n + (a − 1) n , The integral terms in (23)-(29) are G (t) = ∞ 1 g (a) m (a, t) da and G a (t) = ∞ 1 g (a) am (a, t) da. In the case of lipid independence (i.e., g (a) ≡ 1), these integrals evaluate to G (t) = 1 and G a (t) = A M (t) M(t) .

Appendix B: Total System Lipid and Average Lipid Per Cell
The ODE system (25)-(29) can be used to derive additional ODEs for the total amount of lipid in the system L (t), the average lipid content per live cellĀ M (t), and the average lipid content per apoptotic cellĀ P (t). In what follows, we make the definition where F (t) is proportional to the dimensionless macrophage recruitment rate.
For L (t), we have: The mechanisms that add lipid to the system are monocyte recruitment, LDL consumption, and local macrophage proliferation (first, second and third terms on the right-hand side). Contrastingly, the sole mechanism of lipid removal from the system is macrophage emigration (final term on the right-hand side). The time evolution of L (t) is driven entirely by the behaviour of live macrophages, while explicit lipid-dependent effects appear only in the proliferation and emigration terms.
ForĀ M (t), we have: The first three terms on the right-hand side show that LDL consumption, necrotic lipid consumption and efferocytosis all act to increaseĀ M (t). On the other hand, monocyte recruitment and macrophage proliferation (terms four and five) always act to reducē A M (t) (note thatĀ M (t) > 1 in all but the extreme case where m (a, t) is a delta distribution). The final two terms, which relate to apoptosis and emigration, can act to either increase or decreaseĀ M (t) depending on their signs. Letting denote β or γ , we see that these terms act to increaseĀ M (t) when G a (t) G (t) <Ā M (t) and act to decreaseĀ M (t) when G a (t) G (t) >Ā M (t). The ratio G a (t) G (t) can be interpreted as the average lipid load of cells leaving the live cell population via apoptosis or emigration at time t. Thus, when the cells leaving the live cell population have a smaller (larger) average lipid load than the live cells that remain, the average lipid load of the live cell population tends to increase (decrease). The impact of apoptosis and emigration onĀ M (t) observed here is entirely due to the lipid-dependent terms in the model. In the absence of lipid-dependent apoptosis and emigration, the final two terms on the right-hand side of (35) equate to zero.
ForĀ P (t), we have: This equation contains only a term relating to macrophage apoptosis. We see that A P (t) increases when G βa (t) G β (t) >Ā P (t) and decreases when G βa (t) G β (t) <Ā P (t). This makes sense becauseĀ P (t) increases (decreases) when the average lipid load of dying cells is larger (smaller) than the average lipid load of those already dead. Unlike the corresponding term in theĀ M (t) equation, the right-hand side here does not vanish in the absence of lipid-dependent apoptosis. Instead, the term in brackets reduces tō A M (t)−Ā P (t) and the ODE acts to equilibrate the average apoptotic cell lipid content to the average live cell lipid content.

Appendix C: Numerical Solution Methods
The model equations are solved by reformulating (23)-(29) as a large system of coupled ODEs (method of lines) and integrating with the MATLAB routine ode15s. This approach requires the semi-infinite a domain for Eqs. (23) and (24) to be truncated at a finite upper limit a max and appropriately discretised. In choosing a max , we aim to minimise numerical error associated with the loss of live cells (and their ingested