Blood Lipoproteins Shape the Phenotype and Lipid Content of Early Atherosclerotic Lesion Macrophages: A Dual-Structured Mathematical Model

Macrophages in atherosclerotic lesions exhibit a spectrum of behaviours or phenotypes. The phenotypic distribution of monocyte-derived macrophages (MDMs), its correlation with MDM lipid content, and relation to blood lipoprotein densities are not well understood. Of particular interest is the balance between low density lipoproteins (LDL) and high density lipoproteins (HDL), which carry bad and good cholesterol respectively. To address these issues, we have developed a mathematical model for early atherosclerosis in which the MDM population is structured by phenotype and lipid content. The model admits a simpler, closed subsystem whose analysis shows how lesion composition becomes more pathological as the blood density of LDL increases relative to the HDL capacity. We use asymptotic analysis to derive a power-law relationship between MDM phenotype and lipid content at steady-state. This relationship enables us to understand why, for example, lipid-laden MDMs have a more inflammatory phenotype than lipid-poor MDMs when blood LDL lipid density greatly exceeds HDL capacity. We show further that the MDM phenotype distribution always attains a local maximum, while the lipid content distribution may be unimodal, adopt a quasi-uniform profile or decrease monotonically. Pathological lesions exhibit a local maximum in both the phenotype and lipid content MDM distributions, with the maximum at an inflammatory phenotype and near the lipid content capacity respectively. These results illustrate how macrophage heterogeneity arises in early atherosclerosis and provide a framework for future model validation through comparison with single-cell RNA sequencing data.


Introduction
Atherosclerosis is a chronic inflammatory condition of the artery wall (Bäck et al. 2019).The disease begins with the retention of low-density-lipoprotein (LDL) particles in the artery wall.LDL particles, which carry fatty compounds called lipids, enter the artery wall from the bloodstream and are retained via interactions with extracellular matrix.Retained LDL (rLDL) particles are rapidly modified via oxidation and aggregation.The accumulation of rLDL particles triggers an immune response that attracts monocyte-derived macrophages (MDMs) to the lesion.MDMs are typically the most numerous immune cell type in early atherosclerotic lesions (Willemsen and de Winther 2020).They play a key role in disease progression by ingesting extracellular lipid and offloading lipid to high-density lipoprotein (HDL) particles, which also enter the lesion from the bloodstream (Kloc et al. 2020).Importantly, MDMs may adopt a variety of phenotypes depending on their interaction with the lesion microenvironment (Tabas and Bornfeldt 2016;Bäck et al. 2019).This includes inflammatory (M1-like) and resolving (M2-like) phenotypes.Over time, sustained inflammation and the death of lipid-laden MDMs may cause the lesion to transition into an atherosclerotic plaque with a large core of extracellular lipid (Guyton and Klemp 1996;Gonzalez and Trigatti 2017).The rupture of this plaque releases the lipid core into the bloodstream, where it promotes blood clot formation and can induce an acute clinical event.Plaque rupture is the most common cause of myocardial infarction (Costopoulos et al. 2017) and a leading cause of ischaemic strokes (Rothwell 2007).Understanding how the in vivo distribution of MDM phenotype is influenced by MDM lipid content and blood LDL/HDL densities are active areas of research.
Macrophages in atherosclerotic lesions exhibit a continuum of inflammatory to resolving phenotypes (Leitinger and Schulman 2013;Bäck et al. 2019).This view supersedes the traditional dichotomous M1/M2 classification of macrophage phenotype; M1 and M2 polarisation now typically refer to the extremes of a phenotype continuum (Barrett 2020).Macrophage phenotype modulation appears to be reversible (Barrett 2020;Lin et al. 2021;Wang et al. 2014), and is largely determined by the balance between (pro-)inflammatory and (pro-)resolving mediators (Tabas and Bornfeldt 2016;Bäck et al. 2019).Following the classification presented in Tabas and Bornfeldt (2016), inflammatory mediators include cytokines such as TNF and IL-1, that are secreted by MDMs upon uptake of modLDL (Liu et al. 2014), and damage associated molecular patterns (DAMPs) that are released upon the secondary necrosis of apoptotic cells (Sachet et al. 2017).Resolving mediators include the cytokines IL-10 and IL-13, and specialised pro-resolving lipid mediators.Resolving mediators are synthesised by macrophages upon apoptotic cell uptake (Decker et al. 2021) and interaction with HDL (Serhan and Levy 2018).LDL and HDL promote the synthesis of opposing mediator types (inflammatory and resolving respectively) and, so, are likely to induce opposing effects on MDM phenotype.
Of particular relevance is the recent lipid-structured model of Chambers et al. (2023), which serves as the foundation of the present study.By extending this model to account for simultaneous variation in MDM phenotype and lipid content, we provide a mechanistic framework to explore the diversity of lipid-associated macrophage states revealed by single-cell RNA sequencing (Dib et al. 2023).Other authors have proposed dual-structured mathematical models (e.g.Bernard et al. 2003;Doumic 2007;Laroche and Perasso 2016;Hodgkinson et al. 2019, reviewed in Kang et al. (2020)).A key difference between these existing models and ours relates to the time evolution of the structure variables: in most existing models, the structure variables are independent whereas in our model their time evolution is coupled.We use our dual-structured model to address the following questions: 1. How do blood LDL/HDL levels impact lesion composition?a.How do they affect the time-evolution of lesion composition?b.How do they affect lesion composition at steady state? 2. How are phenotype and lipid content distributed among MDMs?a. How do MDM phenotype and lipid content evolve over time?b.Are MDM phenotype and lipid content correlated?c.What are the qualitative features of the phenotype and lipid content marginal distributions at steady state?
The remainder of the paper is structured as follows.Sect. 2 details the model development, including the derivation of a closed subsystem and discussion of parameter values.Sect. 3 contains the results of our model analysis.Key questions 1 and 2 are addressed in Sects.3.1 and 3.2, respectively.Finally, we discuss the results and their implications in Sect. 4.

Model Development
In this section we present a phenotype-lipid dual-structured model for MDM populations in early atherosclerosis.Model schematics for the MDM dynamics and LDL retention are given in Figs. 1 and 2, respectively.

Assumptions and Definitions
We assume for simplicity that macrophage phenotype and lipid content change by finite increments, φ = 1 and a > 0 respectively.Specifically, we let m φ, (t) denote the number density of MDMs with phenotype φ and lipid content a 0 + a at time t ≥ 0. The phenotype index runs over both positive and negative integer values: φ = 0, ±1, . . ., ±φ max .Macrophages with φ > 0 are pro-inflammatory and have an M1-like phenotype; those with φ < 0 are anti-inflammatory and have an M2-like phenotype.The extreme values φ = ±φ max can be interpreted as complete M1 and M2 polarisation.The lipid index runs over non-negative values: = 0, 1, . . ., max , so that macrophage lipid content ranges from their endogenous content, a 0 ≥ 0, to a maximum value, a 0 + κ, where κ := max a is the maximum capacity for ingested lipid.
We also introduce variables to describe the extracellular environment.We denote by L L DL (t) ≥ 0, L r (t) ≥ 0, L ap (t) ≥ 0 and L n (t) ≥ 0 the mass densities of free LDL lipid, retained LDL lipid, apoptotic lipid and necrotic lipid respectively.We let H (t) be the lipid capacity of HDL particles in the lesion.Finally, we denote by S + (t) ≥ 0 and S − (t) ≥ 0 the mass densities of inflammatory and resolving mediators respectively.
We note that all dependent variables in the model denote densities within the tunica intima.This is the innermost region of the artery wall where atherosclerosis develops (Bäck et al. 2019).As shown in Figs. 1 and 2, the intimal densities are nonetheless affected by the bloodstream and tunica media (the next layer of the artery wall) via the flux of cells and lipids.We assume for simplicity that the contents of the bloodstream and tunica media are static in the equations below.), which diffuses freely between the bloodstream, tunica intima and tunica media, and retained LDL (L r ) that is bound to the ECM within the tunica intima.We assume in Eqs. ( 7)-( 8) that the ECM has a finite capacity for rLDL, K r , and that the rate of binding is proportional to the available capacity:

Model Equations
MDMs.We propose that the MDM population evolves according to the following ODEs: apoptosis and egress
The first term on the right hand side of Eq. (1) accounts for lipid uptake.Following Chambers et al. (2023), we model lipid uptake with a mass-action treatment of the following reactions: for every (φ, ) and i ∈ {LDL, r, ap, n}.Reactions (3) assume that the rate of lipid uptake decreases linearly with lipid content in a manner commensurate with available capacity.For notational brevity in Eq. ( 1), we introduce vectors for the uptake rates, k L := (k L DL , k r , k ap , k n ), and the extracellular lipids, The second term of the RHS of Eq. ( 1) accounts for lipid efflux to HDL.Again, following Chambers et al. (2023), we treat efflux with mass-action kinetics according to the reactions: for every (φ, ).The reactions (4) assume that the efflux rate increases linearly with lipid content.The third and fourth terms of the RHS of Eq. ( 1) account for MDM phenotype modulation by inflammatory and resolving mediators (Bäck et al. 2019), S + and S − respectively.We model phenotype modulation via the reactions: for every (φ, ).Here k S is the rate of mediator binding to MDM surface receptors and s is the mediator mass bound per interaction.We assume each mediator binding interaction stimulates a move to a new phenotype class with probability: The parameter 0 ≤ χ ≤ (2φ max ) −1 which modulates this probability can be interpreted as the phenotypic plasticity of the MDMs.We assume that p ± φ decreases linearly to zero as φ → ±φ max , so that it becomes increasingly difficult for MDMs to become more polarised the more polarised they are.Biologically, this property reflects the saturation of intracellular signalling pathways when macrophages are continually exposed to inflammatory or resolving mediators.
The final terms on the RHS of Eq. (1) account for MDM recruitment, apoptosis and egress.We assume that newly recruited MDMs enter the lesion from the bloodstream carrying only endogenous lipid and with an uncommitted phenotype: The recruitment rate is a first-order Hill function of the inflammatory mediator density, S + , which saturates at the maximum value σ M .Resolving mediators, S − , inhibit recruitment by linearly increasing the threshold for half-maximal recruitment from the basal value S c50 + ; the parameter ρ governs the sensitivity of the recruitment rate to S − .We assume for simplicity that the rates of MDM apoptosis, β, and egress, γ , are constant.
Extracellular lipids.We assume that the densities of free and retained LDL evolve according to the ODEs: The first two terms on the RHS of Eq. ( 7) account for the flux of free LDL between the lesion and the lumen/tunica media.The parameters π (0) L and π (1)  L denote the LDL exchange rate at the endothelium and internal elastic lamina, respectively.The densities of free LDL in the lumen, L (0) , and tunica media, L (1) , are assumed to be non-negative constants.We assume that free LDL binds to ECM proteoglycans with rate constant k b in a capacity-limited manner, and unbinds at rate k −b ; the maximum capacity for LDL retention is K r .The final terms on the RHS of Eqs. ( 7) and (8) account for MDM uptake of free and bound LDL via reactions (3).
We propose that the densities of apoptotic and necrotic lipid satisfy: The first term of the RHS of Eq. ( 9) accounts for lipid deposition into the extracellular space due to MDM apoptosis.We assume that apoptotic cells undergo secondary necrosis at rate ν, which provide linear sink and source terms in Eqs. ( 9) and (10), respectively.The final terms on the RHS of Eqs. ( 9) and (10) describe apoptotic and necrotic lipid uptake by MDMs.
HDL lipid capacity.We assume that the lipid capacity of the HDL particles in the lesion evolves according to: The first two terms on the RHS of Eq. ( 11) describe the flux of HDL lipid capacity (via HDL particle diffusion) between the lesion and the lumen/tunica media.For simplicity, we assume that the diffusivity of HDL particles is independent of their lipid capacity, so that π (0) H and π (1) H represent the common permeability of HDL particles and HDL lipid capacity at the endothelium and internal elastic lamina respectively.The final term accounts for MDM lipid efflux to HDL particles.Inflammatory and resolving mediators.We propose that the densities of inflammatory and resolving mediators satisfy the following ODEs: The first term on the RHS of Eq. ( 12) models the release of inflammatory signals by resident cells (e.g.smooth muscle cells, tissue-resident macrophages) (Williams et al. 2019(Williams et al. , 2020)).The signal cascades which first stimulate MDM recruitment are, as yet, unknown, but are thought to be the result of excessive LDL retention (Williams and Tabas 2005).Hence, we assume for simplicity that the resident cells produce inflammatory mediators at rate proportional to the retained LDL lipid density.
Equations ( 12) and (13) also account for mediator production due to lipid activity.This includes MDM production of inflammatory mediators due to LDL uptake (both native and modified forms of LDL induce inflammatory responses in macrophages (Allen et al. 2022;Chen and Khismatullin 2015)), and MDM production of resolving mediators upon apoptotic lipid uptake (Dalli and Serhan 2012) and interaction with HDL (Serhan and Levy 2018).We also include production of inflammatory DAMPs that are released by apoptotic bodies upon secondary necrosis (Sachet et al. 2017).We assume for simplicity that the magnitude of mediator production is proportional to the amount of lipid involved in each of the above interactions; the parameter μ represents the mediator mass produced per unit lipid.
The remaining source terms in Eqs. ( 12) and ( 13) account for constitutive mediator production by MDMs.Indeed, macrophages are potent cytokine emitters (even in the absence of stimulants (Chen and Khismatullin 2015)) and exhibit a phenotypedependent secretion profile (Kadomoto et al. 2021).To account for these effects, we assume that MDMs constitutively produce mediators at the constant rate 2k c per cell, but that the ratio of inflammatory to resolving mediator production is skewed linearly according to phenotype.M1-polarised cells with φ = φ max secrete only inflammatory mediators while M2-polarised cells with φ = −φ max emit only resolving mediators.
We further assume that mediators bind to MDM surface receptors, according to reaction (5), and undergo natural decay at rate δ S .We use a common decay rate for both inflammatory and resolving mediators since experimentally reported half-lives for inflammatory and resolving cytokines are comparable (Liu et al. 2021).Similarly, we use a common binding rate in the absence of evidence for MDM preferential binding.Initial conditions We close Eqs.( 1), ( 7)-( 13) by supposing that at t = 0: The conditions (42) describe the atherosclerotic lesion immediately prior to MDM recruitment.These expressions are derived by solving Eqs. ( 1), ( 7)-( 13) at steady state with m φ, = 0 for every (φ, ).We find that free LDL lipid and HDL lipid capacity are balanced by their fluxes at the endothelium and internal elastic lamina.
Retained LDL levels reflect a balance of binding/unbinding kinetics and directly scale inflammatory mediator levels.The remaining variables are initially zero.
123  123 The blood densities of LDL lipid, L (0) , and HDL lipid capacity, H (0) , are key parameters in our model.Since these quantities are sensitive to modifiable lifestyle factors such as diet and exercise (Schoeneck and Iggman 2021), we explore a range of plausible values in our analysis.The range for L (0) is based on the human serum LDL cholesterol distribution reported in Lee et al. (2012), in which 99.8% of the subjects had LDL cholesterol below 280mg/dL.We multiply this figure by 78/50 to account for LDL phospholipid to obtain an upper estimate of 480mg/dL (Orlova et al. 1999).We take H (0) = 230mg as a conservative upper bound, which exceeds the highest recorded HDL cholesterol concentrations ( 193 mg/dL) in a sample of 116508 individuals from the general population (Madsen et al. 2017).
We consider a range of values for the LDL retention capacity, K r , since it varies between artery wall sections (Lewis et al. 2023).As LDL retention is driven by LDLproteoglycan binding, we estimate K r by considering the artery wall proteoglycan density.Artery wall extracellular matrix prior to atherosclerosis-induced collagen degradation consists of 4% proteoglycan and 40% collagen (Wight 2018).Using 0.01-3 mg/mL as an estimate for the collagen density (0.75 mg/mL is used to replicate the tunica intima in culture models (Liu et al. 2023)), we obtain a proteoglycan density of 0.001-0.3mg/mL.Assuming each proteoglycan molecule (est.mass 800kDa (Yoneda et al. 2002)) can support a single aggregate of LDL (typical diameter 75nm, corresponding to 0.0002 pg lipid for spherical droplets (Guyton and Klemp 1989)), we estimate K r ≈ 15 − 4500 mg/dL.
The values of the remaining parameters are fixed.See Appendix A for further details on these choices of parameter values.

Non-dimensionalisation
We recast the model in terms of the following dimensionless variables: This scaling measures time in units of mean MDM lifespan, β −1 ≈ 1month, and MDM densities relative to the maximum influx per MDM lifespan, σ M β −1 ≈ 17000 mm −3 .We express the extracellular lipid densities and HDL lipid capacity relative to the maximum influx of MDM endogenous lipid per MDM lifespan, a 0 σ M β −1 ≈ 45 mg/dL.Mediator densities are measured relative to the density for half-maximal MDM recruitment, S c50 + ≈ 5 ng/mL.We also introduce a number of dimensionless parameters in Table 2.

Numerical Solutions and Timescales
The dimensionless parameters in Table 2 span several orders of magnitude.In particular, the mediator parameters α, δ S , k S and μ are considerably larger than the other constants.Numerical solutions of Eqs. ( 29)-( 41) consequently require small timesteps to maintain stability.We address this issue in our numerical solutions, computed with Wolfram Mathematica, by using the routine NDSolve with the "StiffnessSwitching" option.Another way to reduce numerical stiffness is to approximate the mediator dynamics via separation of timescales.With δ −1 S 1 and assuming that α, k S , μ = O(δ S ), it is straightforward to show that S ± satisfy the following uniformly-valid quasi-steady state approximations: where α := α/δ S , μ := μ/δ S and kc := k c /δ S .Although approximations ( 43)-( 44) are not used for the simulations presented in Sect.3, they reveal that, at leading order, the mediator densities are proportional to their net rates of production.

Results
We present the results in two sections.In Sect.3.1 we analyse the subsystem (32)-( 41) to generate insight into lesion composition.In Sect.3.2 we focus on the MDM phenotype-lipid distribution, m φ, .

Lesion Composition
We begin our analysis of lesion composition by computing time-dependent solutions of the subsystem (32)-( 41); the results are presented in Sect.3.1.1.In Sect.3.1.2we then analyse the impact of the key parameters: L , H and K r on the steady state values.

Time Evolution
Two typical numerical solutions of the subsystem (32)-( 41) are shown in Fig. 3.The left and right solutions respectively correspond to healthy (L = 3, H = 2.5) and unhealthy (L = 4.5, H = 1) blood levels of LDL lipid and HDL capacity.We set K r = 10 for both cases.The quantity L(t) in plots (d) and (i) is the total extracellular lipid density: 123 whereas plots (e) and (j) show the components of the total lesion lipid content: which includes MDM lipid.The model dynamics can be broadly divided into three phases that we detail below.
The first weeks after initial MDM influx (0 < t < 0.2) are characterised by a decline in lesion LDL and rLDL content (Fig. 3e, j).Both reductions are driven by MDM uptake of rLDL, which also promotes binding of free LDL by restoring the LDL retention capacity.The rise in L M is accompanied by increases in ˆ M and S + since rLDL uptake stimulates production of inflammatory mediators that drive inflammatory phenotype modulation (see Fig. 3b, c, g, h).
For 0.2 < t < 2, the growth of L M and S + slows and ˆ M decreases (see Fig. 3b, c, g,  h).These behaviours arise because the MDM population has ingested enough lipid for the lipid efflux rate to become comparable to that of uptake.For both simulations, the efflux rate exceeds the uptake rate at t ≈ 0.5, where L M attains a local maximum.The corresponding increase in resolving mediators, S − , relative to inflammatory mediators, S + , promotes resolving phenotype modulation.
For 2 < t < ∞, the model tends to a non-zero steady state in a manner sensitive to parameter values.For healthier balances of blood LDL lipid and HDL capacity (L , H ): mean MDM lipid content declines to a small non-zero value and mean MDM phenotype remains negative (i.e.resolving) (Fig. 3b), resolving mediators outbalance inflammatory mediators (Fig. 3c), extracellular lipid levels are low while HDL capacity remains substantial (Fig. 3d), and lesion lipid content decreases to densities below those prior to MDM influx (Fig. 3e).By contrast, when L is sufficiently high relative to H (quantified in Sect.3.1.2):MDM densities are greater (Fig. 3f), mean MDM lipid content increases to a large value and mean MDM phenotype is positive (i.e.inflammatory) (Fig. 3g), inflammatory mediators outbalance resolving mediators (Fig. 3h), HDL capacity is largely exhausted while extracellular lipid levels remain substantial (Fig. 3i), and lesion lipid densities increase to values higher than those prior to initial MDM influx (Fig. 3j).
In Fig. 4 we show how the LDL retention capacity, K r , impacts the timescale of lesion development.The left plot depicts the time to steady state, defined numerically as the smallest time t for which i 1 y i dy i dt 2 ≤ 10 −8 , where the sum is taken over all subsystem variables: y = (M, ˆ M , L M , H , L, S ± ).The time to steady state increases with K r .The trend is more pronounced in the unhealthy case (L = 4.5, H = 1), which exhibits a 20-fold increase over 0.3 ≤ K r ≤ 100.The right plot shows the time, t M , for MDM lipid to exceed 450mg/dL (the smallest t satisfying (1 + κ L M (t))M(t) > 10), which we use as a proxy for fatty streak onset.We find that increases to K r yield smaller values of t M ; regions with higher LDL retention capacity develop fatty streaks earlier.Moreover, if K r is sufficiently small, the MDM lipid density never exceeds the 450mg/dL threshold; fatty streaks will not develop in regions of sufficiently low LDL retention capacity.

Steady State Solutions
The results of Sect.3.1.1show that the model tends to a non-zero equilibrium as t → ∞.Steady state solutions to the subsystem (32)-( 41) are computed numerically via the Mathematica FindRoot routine.
Figure 5 illustrates how blood levels of LDL lipid, L , and HDL capacity, H , impact the model lesion at steady state.Lesion composition becomes more pathological as L increases and H decreases; MDM density, mean phenotype, lipid content and total lipid content each monotonically increase with L and decrease with H , while HDL capacity monotonically decreases with L and increases with H .We note that all plots exhibit regions in which the contours are approximately linear.For the case K r = 10 shown, these contours take the form 0.4L − H = constant.The greater weighting on H in the linear combination reflects the higher value of the dimensionless lipid efflux rate, k H , relative to the uptake rates in Table 2. Values of (L , H ) above the ˆ M = 0 contour (0.4L − H ≈ 0.4) yield healthy lesions with smaller MDM densities, resolving mean phenotypes and small mean lipid loads.Lesion lipid content is low and HDL capacity remains substantial.Values of (L , H ) below the ˆ M = 0 contour yield pathological lesions with higher MDM densities, inflammatory mean phenotypes and higher mean lipid loads.Lesion lipid content is also large.These markers of pathology each increase with 0.4L − H .By contrast the HDL capacity is exhausted in this region.
Figure 6 shows how the LDL retention capacity, K r , impacts lesion composition at steady state.We find that MDM density, MDM mean phenotype and MDM lipid content increase non-linearly with K r .These trends are amplified when the values (L , H ) are more pathological (i.e.0.4L − H is larger).HDL capacity instead decreases with K r .Moreover, total lesion lipid content increases or decreases with K r if the values (L , H ) are pathological or healthy respectively.This difference arises because higher LDL retention capacities give rise to higher MDM densities, due partially to elevated signalling by resident cells.In healthy cases, HDL capacity remains substantial so that MDMs provide a net reduction in lesion lipid content by consuming extracellular lipid and efficiently offloading to HDL.In pathological cases, HDL capacity is exhausted and lipid taken up by MDMs can only leave the lesion via MDM egress.This makes MDMs net contributors to lesion lipid content; on average, they remove less lipid from the lesion than they supply via their endogenous lipid content.

MDM Phenotype-Lipid Distribution
We now turn to Eq. ( 29) to study how phenotype and lipid content are distributed amongst lesion MDMs.After presenting numerical solutions for m φ, (t) in Sec.3.2.1,we derive a continuum approximation of Eq. ( 29) in Sect.3.2.2 to make analytical progress.In Sect.3.2.3we show how analysis of the leading-order advective dynamics enables us to characterise the expected trajectories of MDMs through phenotype-lipid space.Finally, in Sect.3.2.4we analyse the MDM distribution at steady state.

Time Evolution
In Fig. 7 we present the dynamics of the MDM phenotype-lipid distribution, m φ, (t), for the cases shown in Fig. 3).At early times (0 ≤ t ≤ 0.2, the distribution evolves in a wave-like manner, from (φ, ) = (0, 0), in the direction of increasing φ and due to the initial phase of LDL uptake noted in Sect.3.1.1.At later times (0.2 ≤ t ≤ 2), the distribution moves leftwards, becoming concentrated at lower values of φ.This resolving phenotype modulation reflects greater production of resolving mediators, S − , relative to inflammatory mediators, S + , in this phase (c.f.Fig. 3e, h).At longer times, the distribution settles to a steady state.In the healthy case, lipid loads gradually reduce to equilibrium as HDL capacity increases and the extracellular lipid densities decrease to their equilibrium values (c.f.Fig. 3d, e).In the unhealthy case, lipid loads increase towards equilibrium as the extracellular lipid densities increase and HDL capacity declines to their steady state values (c.f.Fig. 3i, j).
We note that the MDM distribution, in both cases and at all times, is concentrated about a central curve which begins at the origin (φ, ) = (0, 0) and terminates at an interior point in (φ, ) space.This feature indicates that phenotype and lipid content are non-linearly correlated.At steady state, this correlation is monotonic but may be negative (as in case a) or positive (as in case b) depending on model parameters.
It is not straightforward to understand how the MDM distributions in Fig. 7 arise by directly considering Eq. ( 29).In order to make progress, we consider a continuum approximation of Eq. ( 29) in the analysis below.
We derive boundary conditions for Eq. ( 49) by requiring that the dynamics of the MDM population, M(t), as defined by Eq. ( 32), are consistent with those of the continuous MDM distribution (48).Integrating Eq. ( 49) over ( φ, ˜ ) ∈ R := [−1, 1] × [0, 1] and applying the divergence theorem yields: where n is the outwards pointing normal vector.We set: Equations ( 53) and ( 54) are no-flux conditions so that MDMs cannot enter or leave the system by exceeding the phenotype bounds φ = ±1 or maximal lipid content ˜ = 1 respectively.We use the Dirac-delta distribution, δ 0 , in Eq. ( 55) to ensure that MDMs enter the lesion at the origin ( φ, ˜ ) = (0, 0).
The dynamics of the velocity field defined by Eqs. ( 57)-( 58) are illustrated in Fig. 8.A striking feature of the plots is that, at all times, the vector field has a single point ( φ× , ˜ × ) of zero instantaneous velocity.This time-dependent point represents the target phenotype and lipid content that MDMs are driven towards by advection.Setting time derivatives to zero in Eqs. ( 57)-( 58), and solving for φ and ˜ shows that: The time-evolution of the zero velocity point, given by Eq. ( 59), is shown in Fig. 8c).We find that ( φ× , ˜ × ) first moves from an inflammatory position towards resolving phenotypes.This initial transition occurs because resolving mediator production by MDMs opposes the initially pure inflammatory environment due to rLDL-stimulated resident cells.When the LDL/HDL balance is healthy, the target phenotype becomes resolving and the target lipid content decreases towards steady state.When the LDL/HDL balance is unhealthy, the target phenotype remains inflammatory and the target lipid content increases as the system evolves to its steady state.We note further that ( φ× , ˜ × ) aligns well with the MDM distribution end points shown in Fig. 7.
Since MDMs enter the model lesion with minimal lipid content and a neutral phenotype, we are particularly interested in solutions to Eqs. ( 57)-( 58) which satisfy the initial conditions: Such solutions represent the trajectory of MDMs that enter the lesion at time t = t in the limit → 0. They can be expressed in terms of the subsystem variables as follows: where I 1 (t) and I 2 (t) are integrating factors: 123 We use numerical solutions of the subsystem ( 32) -( 44) to evaluate Eqs. ( 61) -( 62) in Fig. 9.As MDMs travel along the trajectories shown, they emigrate from the lesion and die via apoptosis at the constant rate (1 + γ ).The probability that MDMs which enter the lesion at time t = t are still alive and in the lesion at time t ≥ t is given by e −(1+γ ) (t−t ) .This probability is represented in Fig. 9 by the opacity of the curves.The plots illustrate that MDMs which enter the lesion early (e.g. at t = 0.2) can be expected to first transition from phenotypic neutrality to an inflammatory state with moderate lipid loads, before evolving to either a resolving phenotype with low lipid load (case a), or a milder inflammatory phenotype with high lipid load (case b).The trajectories of MDMs that enter the lesion at later times are typically monotonic, in contrast to the looping trajectories of MDMs that enter at earlier times.We note in particular that trajectories with t = 100 (near steady state) align well with the centre line of the MDM distributions (c.f.Fig. 7 at t = 100).

Steady State Solutions
It is straightforward to show that m ( φ, ˜ ), the steady state solution of Eq. ( 49), satisfies: 123 where the constants φ∞ , ˜ ∞ , q and p are defined in terms of model parameters and the equilibrium values of the subsystem variables: We note that φ∞ ∈ (−1, 1) and ˜ ∞ ∈ (0, 1) coincide exactly with the target phenotype and lipid content, φ× and ˜ × respectively (c.f.Eq. ( 59) at steady state).We note also that p ∈ (−1, ∞) and q ∈ (0, ∞).The boundary conditions ( 53)-( 55) become: where the dimensionless fluxes are given by: Rather than searching for a closed-form solutions to Eq. ( 63) in full generality, we characterise the solutions via asymptotic analysis in the limit → 0.More specifically, we derive equations for the central curve of the MDM distribution and marginal distributions with respect to lipid content, W ( ˜ ), and phenotype, V ( φ), where: Finally, we consider the impact of the LDL lipid and HDL capacity blood densities, L and H respectively.The results are summarised in Figs. 10 and 11.
Central curve We determine the curve about which m ( φ, ˜ ) is centred by considering Eq. ( 63) at leading order.Indeed, when = 0, Eq. ( 63) reduces to a first order hyperbolic PDE, with characteristics that satisfy the equation: Solving Eq. ( 69) subject to φ( ˜ = 0) = 0 yields: (same colour legend as in Fig. 7).The second and third columns compare the exact solutions for the lipid content and phenotype marginals, given by Eqs. ( 73) and ( 83) respectively, with the discrete marginals.We set K r = 10 for all cases (Color figure online) In the first column of Fig. 10, the solution (70) is superimposed on numerical solutions for m φ, at steady state, showing good agreement with the centre of the MDM distributions.Equation (70) confirms that, at steady state, phenotype and lipid content are monotonically correlated via a power law.The correlation is positive when φ∞ > 0 and negative when φ∞ < 0. Equivalently, using the first of Eqs. ( 64), if inflammatory mediators dominate resolving ones (S + > S − ) then MDMs with higher lipid loads have a more inflammatory phenotype; if, however, resolving mediators dominate (S + > S − ) then MDMs with higher lipid loads have a more resolving phenotype.The nonlinearity of the correlation is determined by the constant q, which measures the relative amount of mediator-MDM activity to lipid-MDM activity.Distribution of lipid content Integrating Eqs. ( 63) and ( 65) with respect to φ ∈ [−1, 1] yields the following ODE for W ( ˜ ): and the boundary conditions: Equation ( 71) admits an exact solution (found via the Mathematica DSolve routine): where f is the function: U b a (z) is the confluent hypergeometric function and L b a (z) the generalised Laguerre polynomial.The constants K 1 and K 2 are determined by substituting Eq. ( 73) into the boundary conditions (72); exact expressions are readily obtained via the Mathematica Solve routine, but are too involved to be insightful and so are omitted here for brevity.We compare the solution (73) to numerical solutions of the discrete model in the second column of Fig. 10, showing excellent agreement.
The form of Eq. ( 73) makes it difficult to understand how W ( ˜ ) depends on the constants defined in Eq. ( 64).Nonetheless, analysis of the leading order "outer" solution shows how the qualitative behaviour of W ( ˜ ) depends on p and ˜ ∞ .Here we follow the asymptotic analysis in Chambers et al. (2023).Setting = 0 in Eq. ( 71) admits the following solution: where K 3 ≥ 0 is a constant of integration.Equation ( 75) admits three possible behaviours as ˜ → ˜ ∞ from below.For p > 0, W outer decreases monotonically to zero (with non-zero derivative if p < 2); for p = 0, W outer is constant; and for p < 0, W outer diverges to +∞.The full solution does not exhibit such discontinuities at ˜ = ˜ ∞ due to the regularising effects of the second derivative term.The corresponding behaviour for W ( ˜ ) is as follows.For p > 0, W decreases monotonically and smoothly to zero near ˜ = ˜ ∞ (Fig. 10b); for p = 0, W takes a quasi-uniform sigmoidal profile with a rapid decrease to zero near ˜ = ˜ ∞ (Fig. 10c); and for p < 0, W increases monotonically before attaining a local maximum near ˜ = ˜ ∞ (Fig. 10a,  d).In the above, "near ˜ = ˜ ∞ " means an O( 12 ) neighbourhood of ˜ = ˜ ∞ .This scaling can be derived by searching for an inner variable ˆ = ( ˜ − ˜ ∞ )/ n for some ˆ = O(1) in Eq. ( 71); an exponent n = 1 2 is required to bring the second derivative into the dominant balance.Distribution of phenotype Integrating Eqs. ( 63) and ( 65) with respect to ˜ ∈ [0, 1] gives the following boundary value problem for V ( φ): where r := −1 + p+1 q .As in the derivation of a Green's function, we first recast the Dirac-delta source, which describes MDM recruitment, as a jump condition.That is, we seek a solution of the form: where V ± solve the ODE: and satisfy: The exact solution for Eq. ( 79) can be written as: where: Expressions for the constants J ±,1 , J ±,2 can be found by substituting Eq. ( 83) into conditions ( 80)-( 82).As for W ( ˜ ), they are too complicated to be insightful and are omitted for brevity.The comparison of solutions ( 83) with the discrete model output in column 3 of Fig. 10 shows good agreement.We note that the jump condition amounts to a small reduction in slope at φ = 0 that is almost imperceptible for the parameter values in Table 1.
Following the analysis of W ( ˜ ), we compute the leading order "outer" solution to Eq. ( 76): Equation ( 85) is of the same form as Eq. ( 75); φ∞ and r play the same roles for V as ˜ ∞ and p play for W .The main difference, as shown in Fig. 10, is that only the case r < 0 manifests in numerical solutions.Consequently, V ( φ) always attains a local maximum in a small (again O( 12 )) neighbourhood of φ = φ∞ .This observation is supported by the numerical results below.Impact of blood LDL lipid and HDL capacity levels.We conclude our steady state analysis by examining how the qualitative form of the MDM distribution changes as L and H vary.The results are summarised in Fig. 11 where we plot the parameter groupings in Eqs.(64) against L and H ; the four cases shown in Fig. 10 are also indicated.Recall that p and r determine the overall shape (e.g.whether they attain a local maximum) of the lipid and phenotype marginal distributions, W ( ˜ ) and V ( φ) respectively, while φ∞ , ˜ ∞ and q determine the central curve.
We note first how r and p vary with L and H .The value of r decreases monotonically with L and is less sensitive to H . Importantly, we find that r < 0 for all L ∈ (0, 10) and H ∈ (0, 5), indicating that the phenotype marginal attains a local maximum near φ = φ∞ across the range of physiologically plausible values for blood LDL lipid and HDL capacity.By contrast, p ≥ 0 for values of (L , H ) sufficiently close to the origin.This suggests that monotone decreasing and sigmoidal lipid marginals only occur when blood levels of LDL lipid and HDL capacity are both sufficiently low; otherwise the lipid marginal attains a local maximum near ˜ = ˜ ∞ .
The target phenotype, φ∞ , and lipid content, ˜ ∞ , both increase with L and decrease with H .In particular, their plots exhibit the linear contours prominent in Fig. 5; increases to the relative level of LDL lipid to HDL capacity raise both these markers of pathology.
Finally, we note that the exponent q of the central curve is largest when L and H are comparable.This causes the central curve to appear straighter (i.e.φc ( ˜ ) ≈ φ∞ except for the smallest values of ˜ ) in cases of intermediate pathology.Indeed, in

Discussion
In this paper we have developed a mathematical model for early atherosclerosis in which the MDM population is structured by phenotype, φ, and lipid content, .This framework allows for incremental changes in phenotype and lipid content, which contrasts their treatment as binary variables in much of the existing modelling literature.The model couples the MDM dynamics to the densities of: free LDL lipid L L DL (t), retained LDL lipid L r (t), apoptotic lipid, L ap (t), necrotic lipid, L n (t), HDL capacity, H (t), inflammatory mediators, S + (t), and resolving mediators, S − (t).These variables form a closed subsystem of ODEs when coupled with the MDM density, M(t), mean phenotype, ˆ M (t), and mean lipid content, L M (t).This subsystem can be solved independently of the structured MDM population, m φ, (t).
We parameterised the model using data from the biological literature.Where possible, we used human in vivo data (e.g.blood measurements for L (0) and H (0) ) or ex vivo data (e.g.surgical data for π (0) L and postmortem data for L (1) and H (1) ).However, the majority of the model parameters are calibrated to in vitro experiments.We prioritised studies with human cell lines (e.g. for a 0 , k L DL , k r , k H , S c50 + , ρ), and used data from nonhuman cell lines when necessary (e.g.murine data for κ, k ap and k n ).Since the point estimates in Table 1 are likely to carry high degrees of uncertainty,

123
we have cautiously interpreted the results of the current study by focusing on trends rather than precise quantitative outputs.
Our model analysis focused on the impact of three dimensionless parameters: L , H and K r .The quantities L and H are proportional to the blood densities of LDL lipid and HDL capacity, respectively.These vary considerably according to genetic and lifestyle factors.We explored values in the range L ∈ (0, 10) and H ∈ (0, 5); the upper bounds correspond to densities of 450 and 225mg/dL, respectively.The quantity K r is proportional to the capacity for LDL lipid retention.This varies according to the specific region of artery wall under consideration.We explored values K r ∈ (0.3, 100), which span athero-resistant regions with capacity 15mg/dL to athero-prone regions with capacity 7500 mg/dL.
We discuss our findings below in relation to the key questions posed in Sect. 1. Q1.a.How do blood LDL/HDL levels affect the time-evolution of lesion composition?Time-dependent numerical solutions revealed that the model lesion evolves in three phases.The first phase consists of the initial influx of MDMs and corresponding decline in lesion LDL and rLDL content.MDM lipid loads are small but increasing and MDM mean phenotype is inflammatory.The second phase is characterised by modulation towards resolving MDM phenotypes and a slower rise of MDM lipid content.These phases are consistent with observations of macrophage behaviour during acute inflammation.Macrophages first adopt inflammatory phenotypes and transition to resolving phenotypes during a phase of inflammatory resolution and tissue repair (Pérez and Rius-Pérez 2022).To the best of our knowledge, the present model is the first to capture this transition as an emergent property of the dynamics.Rather than completely resolving, however, the model lesion enters a final phase in which the dynamics tend to a nonzero steady state.If blood LDL lipid is low relative to HDL capacity, the lesion tends to a healthy state with low lipid burden and resolving phenotypes.This behaviour suggests that the well-documented spontaneous regression of many early atherosclerotic lesions (Insull 2009) may simply be the natural and expected progression under a healthy blood lipoprotein balance.By contrast, when blood LDL lipid is high relative to HDL capacity, MDM lipid loads increase to equilibrium and the model lesion accumulates necrotic lipid.Overall, the three-phase dynamics support the idea that chronic inflammation in atherosclerosis can be understood as an acute inflammatory response to LDL retention with incomplete resolution (Sansbury and Spite 2016).
We also studied the impact of LDL retention capacity on the timescale of atherosclerosis development.In particular, we computed the time for MDM lipid to exceed a density of 450mg/dL; we use this lipid density as a proxy for fatty streak formation, which is characterised by the appearance of foam cells (Daskalopoulos et al. 2015).Our results indicate that this time decreases with LDL retention capacity; fatty streak onset in the model lesion occurs earlier for regions of high LDL retention.This finding is consistent with observations of the murine aortic arch in which regions of lower retention capacity (e.g.central zone of the arch) developed atherosclerosis later than regions of high retention (e.g.dorsal and ventral zones) (Lewis et al. 2023).

Q1.b. How do blood LDL/HDL levels affect lesion composition at steady state?
Analysis of steady state solutions revealed how equilibrium lesion composition depends on the parameters L , H and K r .The results indicate that the degree of pathology is largely determined by a linear combination of the form: 0.4L − H ; it is the (weighted) relative value of LDL lipid to HDL capacity in the blood that matters.The greater weighting on H in the linear combination reflects the higher value of the dimensionless lipid efflux rate, k H , relative to the uptake rates in Table 2. Hence, the increased amount of lipid efflux promoted by a rise in blood HDL lipid capacity is greater than the increase in lesion lipid content by an equal rise in blood LDL lipid density.Overall, the model predicts that early atherosclerotic lesions regress upon blood LDL lipid density increases or blood HDL capacity increases.The degree of pathology also generally increases with K r .Regions of higher LDL retention capacity exhibit greater MDM densities, more inflammatory MDM phenotypes and higher MDM lipid loads.Q2.a.How do MDM phenotype and lipid content evolve over time?We analysed the MDM distribution by deriving a continuum analogue of the discrete equations (29) for m φ, (t).Analysis of the continuum model showed that the time-evolution of phenotype and lipid content for individual MDMs depends on the time of entry into the lesion.MDMs which enter the lesion at early times are expected to first transition from phenotypic neutrality to an inflammatory state, and then to a resolving state (if the blood LDL-HDL balance is healthy) or to a milder inflammatory state (if the LDL-HDL balance is unhealthy).The trajectories of MDMs which enter the lesion at later times are monotonic (in phenotype-lipid structure space) and follow the central curve of the equilibrium MDM distribution.These results suggest that MDMs which enter the lesion during early stages of lesion development experience a greater amount of phenotype modulation throughout their lifespan than those entering at later times.Q2.b. Are MDM phenotype and lipid content correlated?The asymptotic analysis presented in Sect.3.2.4showed that MDM phenotype and lipid content are correlated via a power law at steady state.If LDL lipid density dominates HDL capacity in the blood, lipid-laden MDMs have a more inflammatory phenotype than lipid-poor MDMs, while if blood LDL lipid density is sufficiently low, lipid-laden MDMs have a more resolving phenotype than MDMs with a lower lipid burden.The non-linearity of the correlation is determined by a constant, q, that measures the relative amount of lipid activity to mediator activity in the lesion (made precise by Eq. ( 64)).Although we did not pursue a time-dependent mathematical analysis, numerical solutions (e.g. Figure 7) demonstrated that the MDM phenotype-lipid distribution is always concentrated about a central curve, indicating that a monotone correlation between phenotype and lipid content holds for all times.These findings are consistent with the recent discovery of PLIN2 hi /TREM1 hi macrophages in human lesions, for which the transcriptional signatures of lipid loading and inflammation are coupled (Dib et al. 2023).Q2.c.What are the qualitative features of the phenotype and lipid content marginal distributions at steady state?Further analysis in Sect.3.2.4showed that the phenotype marginal distribution always attains a single local maximum.The location of the maximum, which represents the most common phenotype in the lesion, is a close approximation to the "target" phenotype that MDMs are driven towards by the extracellular environment over their lifetime.By contrast, the lipid marginal distribution varies more in shape as L and H are varied; it may exhibit a local maximum, adopt a quasi-uniform profile, or decrease monotonically according to the value of p, a constant which quantifies the amount of MDM-lipid activity in the lesion.a typical blood monocyte density (Nelson 2014), gives 34000 mm −3 per month.Since approximately half the incoming monocytes do not differentiate into macrophages and leave the lesion within 48 hours (Lee et al. 2019), we estimate σ M = 17000 mm −3 per month.
The uptake and efflux rates are obtained from in vitro experiments.We note that the free LDL uptake rate, k L DL , is calibrated to data for native LDL (data for oxidised LDL yields a slightly higher but comparable estimate), while the larger rLDL uptake rate, k r , is based on aggregated LDL data.This is reasonable since proteoglycan-bound LDL aggregates into lipid droplets that are more readily taken up by MDMs in vivo (Öörni and Kovanen 2021).The apoptotic lipid uptake rate, k ap , has been scaled down from its in vitro estimate by a factor 19 to account for the 19-fold deficiency in efferocytosis measured in human atherosclerotic plaques (Schrijvers et al. 2005).
The mediator mass, s, and natural decay rate, δ S , are based on data for the inflammatory mediator TNF, which has a molecular weight of 17kDa in its secreted form as a monomer and a 19 minute half-life (Atzeni and Sarzi-Puttini 2013;Liu et al. 2021).Other mediators have a comparable mass and half-life (e.g. 15 kDa and 20 min respectively for the resolving mediator IL-4 (Liu et al. 2021)).
We estimate the rate of rLDL-stimulated mediator production from resident cells, α, using reports that MDMs are observed in murine atherosclerotic lesions within 2 weeks of high fat diet treatment (Williams et al. 2020).The value α = 3.5 × 10 −3 per month ensures that the MDM density in our model is approximately 350 mm −3 (approximately blood monocyte concentration) in 2 weeks when L (0) = 750 mg/dL.
Finally, we note that the lipid uptake/efflux increment, a, and maximal structure indices, max and φ max , are model abstractions.In reality, lipid increments vary substantially depending on the specific interaction.The smallest increment is likely that associated with efflux, given that each HDL particle can store approximately 75kDa of cholesterol.The largest increment is likely that due to apoptotic or necrotic lipid uptake, which we can assume to be smaller than the endogenous lipid content, a 0 , since cellular uptake in vivo is typically piecemeal rather than involving whole cell engulfment (Taefehshokr et al. 2021).We fix a = 16pg as a reasonable intermediate value.This choice means that max = κ/ a ≈ 100.The value of φ max has an upper bound of (2χ) −1 ≈ 83000 since the phenotype transition probabilities in Eq. ( 5) must satisfy p ± φ ≤ 1.We fix φ max = 50 so that phenotype and lipid content are equally resolved in our model.

Fig. 1
Fig. 1 A schematic of the MDM-lipid dynamics in the model.Processes represented with blue arrows stimulate the secretion of resolving mediators by MDMs, while those represented by red arrow stimulate the emission of inflammatory mediators.The lower half illustrates the discrete phenotype-lipid structure space that underpins the MDM dynamics.LDL retention and constitutive mediator production by MDMs are not shown Tarique et al. (2015) a MDM uptake/efflux increment of lipid 16 pg Kontush et al. (2007); Taefehshokr et al. (2021), ∈ (75kDa, a 0 ) max Maximum MDM lipid capacity per uptake/efflux increment 100 = κ/ a φ max Half the maximum number of MDM phenotype classes 50 < (2χ) −1

Fig. 3 123 Fig. 4 Fig. 5
Fig. 3 Time evolution of lesion composition.The plots a-j show numerical solutions of the subsystem (32)-(41) for a case with healthy LDL-HDL balance: L = 3, H = 2.5 (left), and unhealthy LDL-HDL balance: L = 4.5, H = 1 (right).The system tends to a non-zero steady state as t → ∞ with values that depend sensitively on L and H .We set K r = 10 for both cases (Color figure online)

Fig. 7
Fig. 7 MDM phenotype-lipid distribution dynamics for two typical simulations.The plots show numerical solutions of Eq. (29) for m φ, (t) for a case with healthy LDL/HDL balance.a L = 3, = 2.5, and unhealthy LDL-HDL balance.b L = 4.5, H = 1.We use K r = 10 for both simulations.At steady state (see t = 100), the distribution may skew towards resolving (case a) or inflammatory (case b) phenotypes (Color figure online)

Fig. 8 Fig. 9
Fig. 8 Dynamics of the MDM phenotype-lipid velocity field.The vector field plots illustrate the right hand side of Eqs.(57)-(58) using the numerical solutions of the subsystem (32)-(41) shown in Fig. 3.The velocity magnitude is indicated by larger arrows and brighter colours.The point of instantaneous zero velocity is indicated by a circle at each time point in the vector field plots, and plotted against time in the bottom plot.We fix K r = 10 and use L = 4.5, H = 1.0 for case a) and L = 3.0, H = 2.5 for case b) (Color figure online)

Fig. 10
Fig. 10 Comparison of continuum analytical results to numerical solutions for the discrete MDM distribution, m φ, , at steady state.Cases a-d use the values (L , H ) = (3, 2.5), (1.7, 0.8), (1.9, 0.6) and (4.5, 1) and are ordered roughly by pathology.The first column overlays Eq. (70) for the central curve with m φ, (same colour legend as in Fig.7).The second and third columns compare the exact solutions for the lipid content and phenotype marginals, given by Eqs.(73) and (83) respectively, with the discrete marginals.We set K r = 10 for all cases (Color figure online)