The LDL-HDL Profile Determines the Risk of Atherosclerosis: A Mathematical Model

1Atherosclerosis, hardening of the arteries, is the leading cause of death in the United States, and worldwide. The disease triggers heart attack or stroke, with total annual death of 900,000 in the United States and 13 million worldwide.

2Atherosclerosis is a disease in which a plaque builds up inside the arteries. A plaque contains low density lipoprotein (LDL), macrophages, smooth muscle cells (SMCs), platelets, and debris. The plaque constricts the lumen of the blood vessel thereby increasing the shear force of blood flow,. As the plaque continues to grow, the increased shear force may cause rupture of the plaque, possibly resulting in the formation of thrombus (blood clot),, ischemic stroke, and heart attack –.

3The process of plaque development begins with a lesion in the endothelial layer, allowing LDL, to move from the blood into the intima and becoming oxidized LDL (ox-LDL) by free radicals (FRs). FRs are oxidative agents continuously released by bio-chemical reactions within the body, including the intima –. Endothelial cells, sensing the presence of ox-LDL, secrete monocyte chemoattractant protein (MPC-1),, which triggers recruitment of monocytes into the intima. After entering the intima, monocytes differentiate into macrophages, which have an affinity for the ox-LDL –. The ingestion of large amounts of ox-LDL transforms the fatty macrophages into foam cells,. Foam cells secrete chemokines which attract more macrophages,,. SMCs from the media move into the intima by chemotactic forces due to MCP-1,, and platelet-derived growth factor (PDGF),, as well as by haptotaxis by the extracellular matrix (ECM). PDGF is secreted by macrophages, foam cells and SMCs,. ECM is remodeled by matrix metalloproteinase (MMP) produced by a variety of cell types including SMCs, and is inhibited by tissue inhibitor of metalloproteinase (TIMP) produced by macrophages and SMCs. Interleukin IL-12, secreted by macrophages and foam cells,,, contribute to the growth of a plaque by activating T cells,,. Indeed, the activated T cells secrete interferon IFN-γ, which in turn activates macrophage in the intima,,. At the same time that LDL enters the intima, high density lipoprotein (HDL) also enters into the intima, and becomes oxidized by free radicals,. However, oxidized HDL (ox-HDL) is not ingested by macrophages. HDL helps prevent atherosclerosis by removing cholesterol from foam cells, and by the limiting inflammatory processes that underline atherosclerosis. Furthermore, HDL takes up free radicals that are otherwise available to LDL. Some of the key players in the atherosclerosis process are shown in Fig. 1.

4It has long been recognized that the cholesterol concentrations in the blood are indicators of the probability that a plaque will develop: higher LDL and lower HDL concentrations indicate a higher probability of plaque development. Public health guidelines in the U.S. specify what levels of LDL are low risk and what levels are high risk; they also specify what levels of HDL are poor and what levels are near ideal,. However, what is more relevant is to specify the risk associated with combined levels of LDL and HDL, and this is what the present paper addresses. A schematic of the network of atherosclerosis is given in Fig. 2. In this paper, we developed a mathematical model of plaque formation by a system of partial differential equations based on Fig. 2. The aim of the model is to determine the risk of plaque formation for combined levels of LDL and HDL. In particular, we created a “risk-map” for plaque development in the LDL-HDL coordinate plane, where the first quadrant of the plane was divided into regions of high risk, low risk and no risk. Anti-cholesterol drugs are aimed at lowering high levels of LDL, but some drugs are known to also increase the level of HDL. Hence such a risk-map may be important when evaluating the extend to which an anti-cholesterol drug can reduce the risk of atherosclerosis for particular individuals.

5In this paper, we present a mathematical model based on the network shown in Fig. 2. The model includes the variables listed in Table 1. We assume that all cells are moving with a common velocity u; the velocity is the result of movement of macrophages, T cells and SMCs into the intima. We also assume that all species are diffusing with appropriate diffusion coefficients. The equation for each species of cells X has a formwhere the expression on the left-hand side includes advection and diffusion, and FX accounts for various growth factors, bio-chemical reactions, chemotaxis and haptotaxis. The equation for the chemical species are the same but without the advection term. Fig. 3 shows a 2D cross section of a blood vessel with plaque Ω, and a planar cross section of a plaque in the direction along a blood vessel.

6The distribution of LDL, HDL, ox-LDL and free radicals in the intima are described using reaction-diffusion equations, where kL and kH are reaction rates of oxidization, and is the reduction rate of ox-LDL due to ingestion by macrophages. Eqs. (1) and (2) model the evolution of LDL and HDL concentrations. It is assumed that LDL and HDL are lost by reaction of oxidation with free radicals. Equation (3) models the production of ox-LDL due to LDL oxidation by reaction with the radicals (first term on right-hand side) and a reduction of ox-LDL through ingestion by macrophages (second term on right-hand side). Equation (4) models the evolution of free radicals concentration with baseline growth r 0.

7The evolution of macrophage density is modeled byHere the first term on right-hand side accounts for recruitment of macrophages by MCP-1, and the second term accounts for the activation of macrophages by IFN-γ,,.

8The MCP-1 equation is given bywhere the first term on the right-hand side is the production of MCP-1 by endothelial cells, whose density is assumed to be constant, under the influence of ox-LDL.

9The density of T cells, which are primarily CD4+ T cells, satisfies the equationIn this equation, we assume that T cells are activated by IL-12 in conjunction with MHC-II (major histocompatibility complex, class II). Actually, T cells are also activated by IL-1 and IL-6 produced by macrophages and SMCs,,. However, because of lack of experimental data, we do not include the IL-1 and IL-6 explicitly but instead consider their effect implicitly in estimating the parameter. For simplicity, we include the anti-inflammatory effect of IL-10 produced by macrophages only implicitly, by the factor.

10The dynamics of IFN-γ concentration is modeled bywhere the first term on right-hand side represents production of Iγ by T cells.

11The equation of the SMCs density is given byThe first two terms on right-hand side account for chemotaxis by MCP-1,, and PDGF,, and the last term accounts for haptotaxis by ECM.

12The concentration of IL-12 is modeled byThe first term of right-hand side is the production of I 12 by macrophages enhanced by Iγ and resisted by HDL. The production of I 12 by macrophages is resisted by I 10 (which, for simplicity, is accounted by the factor ),. The second term represents the production of I 12 by foam cells.

13We have the following sets of reaction diffusion equations for the chemokines (G, Q and Qr): In equation (11), PDGF is produced by macrophages, foam cells, and SMCs,. In equation (12), MMP is secreted by SMCs (first term on right-hand side), and is lost by binding with TIMP (second term). In equation (13), TIMP is produced by SMCs and macrophages.

14Macrophages that have ingested a large amount of ox-LDL become foam cells,,, so we have

15We assume that the intima has the constituency of a porous medium. Then, by Darcy's law, the velocity u of the cells is given bywhere σ is the pressure. We also assume that the total density of all the cells plus the concentration of ρ is constant. This constant should be smaller than the average density of a plaque, 1.22±0.03 g/cm 3, because plaques contain some debris, which are not included in our model. We take the constant to be 1 g/cm 3, i.e.,We assume that all cells are approximately of the same volume and surface area, so that the diffusion coefficients of the all cells have the same coefficient, D. By adding Eqs. (5), (7), (9) and (14), we getwhere Eq. (17) gives a relation between ρ and σ. We next derive an equation for ρ. The ECM is degraded by MMP, and is remodeled by macrophages and SMCs,. For simplicity, we take the remodeling rate to be a constant, λρ, as in. Then the equation of the density of ECM is given bySince this equation can be written in the formwhere

16Adding this equation and (17), we get our final equation for σ:The equation for can then be written as

17For simplicity, we consider only 2-dimensional plaques as in Fig. 3. Then the boundary of the plaque consists of i) ΓM, in contact with media; ii) a free boundary ΓI, inside the lumen, and iii) two more vertical boundaries ΓL and ΓR of the intima in the case of Fig. 3(B).

18We assume flux boundary conditions of the formfor X = L, H, M, T, and non-flux boundary conditions for all other variables,where, r, P, Iγ, S, I, I 12, G, Q, Qr, F. The boundary values for ρ are determined by Eq. (16). The coefficient is a constant except for M, and, since ox-LDL attracts monocytes, while HDL limits the inflammation process. Note that L 0 and H 0 are the LDL and HDL concentrations in the blood, so we shall be interested to see how these concentrations determine whether a small plaque will grow or shrink.

19As in –, we assume that the free boundary ΓI is held together by cell-to-cell adhesion forces so thatwhere κ is the mean curvature of the surface ΓI. (If ΓI is circular, then κ is the reciprocal of the radius) Furthermore, the continuity condition, where n is the outward normal and Vn is the velocity of the free boundary ΓI in the direction n, yields the relation

20We assume non-flux boundary conditions for all variables except ρ and S on ΓM: For S, we havewhere S 0 is SMCs density in the media, and for simplicity we take since MCP-1 and PDGF attract SMCs from the media. As in the case of ΓI, the boundary values of ρ are determined by Eq. (16).

21We assume the periodic boundary conditions on ΓL and ΓR.

22Table 2 lists the range of molecular weights of proteins and Table 3 lists their range of concentration. In the second columns in Tables 2 and 3, we indicate the (intermediate) values used in the simulations. The Tables 2 and 3 are used to estimate some of the model parameters. A summary of all the model parameters is given in Tables 4 and 5.

23To estimate some of the parameters in the equations for proteins, we shall use the concept of “accessible surface area”, of a protein p, or briefly “area,”, which is roughly the minimum surface area of the smooth shapes containing the protein. It was estimated in, that, and, so that their ratio isAccordingly, the corresponding reaction rates of the oxidation, kL and kH, are related by. Moreover, the reaction rate of oxidation of LDL by free radicals is day−1,,, so that g cm −3 day−1.

24We assume that the diffusion coefficients of all the cells are the same, and take them to be cm 2 day−1,. In order to estimate the diffusion coefficients of the various proteins, we assume that the diffusion coefficient of protein p, Dp, is proportional to its area Ap, i.e.,, where we take K to be the same for all small molecules. For glucose, which is a monomeric globular protein, can be computed in terms of the molecular weight, by the formula, and dalton, day−1. Hence, for glucose, K is determined byWe can now compute day−1 and day−1. Free radicals are monomeric globular proteins (average weight is 500 da, Table 2). HenceThe diffusion coefficient of MMP is day−1. We assume that the diffusion coefficient of TIMP is same as that of MMP.

25We assume that, in Eq. (7),, where denotes the average concentration of X. We take gcm−3, gcm−3 from Table 3, and day−1,. Then is estimated by day−1.

26PDGF is produced by SMCs, and likely also by endothelial cells and macrophages. In wound healing, macrophages produce PDGF at rate of 5.76 day−1. Since the plaque formation is a much slower process, we take this rate to be much smaller, i.e., day−1. Since SMCs production rate of PDGF is higher than that by macrophages, we take day−1.

27The production rate of MMP by tumor cells was estimated in to be day−1. We assume that SMCs produce MMP at a much lower rate, namely, day−1. Since SMCs produce MMP to enable them move into the intima by haptotaxis, we assume that they produce TIMP at a lower rate than MMP, and take day−1. As macrophages produce most of the TIMP, we take the production rate of TIMP by macrophages to be day−1.

28We assume that the production rate of MCP-1 by endothelial cells,, is twice that of, where is the concentration of MCP in the blood, which is equal to. We assume that day−1, and that, in Eq. (14), and, so that day−1.

29By, day−1. We assume that foam cells have lower production rates of I 12 and PDGF than macrophages, and take and to be one third of the values of and, respectively, so that day−1, and day−1.

30The degradation rate of MMP is day−1. Since TIMP has a short half life compared to MMP, we take its degradation rate to be day−1.

31In, the binding rate of MMP and TIMP is reported to be, where the mass per mole, and the molecular weights of MMP and TIMP are 52 kda and 25 kda, respectively. Accordingly, we derive the binding rate per Molar per second (by same formula as in ),andwhere NA is called the Avogadro number, and is the number of molecular per dm 3., and is the mass of a proton for atomic mass unit.

32The range of macrophages in the blood is; we take. The range of T cells in the blood is; we take. The range of SMCs is; we take. We assume that is half of in Eq. (6), and similarly, in Eq. (7), and in Eq. (10). We assume that the influx of LDL and HDL into the intima is larger than the influx of macrophages and SMCs, and take, and. The influx of T cells is assumed to be smaller than that of macrophages, and we take.

33In order to illustrate our numerical method, we consider the following diffusion equation with Robin boundary conditions: where is an advection term, and either or (no advection), and. Multiplying the differential equation by an arbitrary function, and performing integration by parts using the boundary conditions, we get

34This is an equivalent formulation of the system (27), which is better suited for simulation.

35Similarly, Eq. (21) for σ has the equivalent form:for an arbitrary function.

36The standard Galerkin discretization method uses finite dimensional subspaces to approximate the solution X. Let be a basis of, where N is the number of nodes within the triangulation K. Let denote the numerical approximation of X at time, where dt is the time step, is written as a sumfor coefficient to be determined. If is approximated by, then (28) is equivalent toor,

37Recalling (30), we can rewrite the system (32) as a linear system of equationswhere is the vector of, and the coefficient matrix and the right-hand side are defined byandSimilarly, setting, where is a numerical approximation of σ at time ndt, Eq. (29) can be written as followswhere is the vector, B is a matrix,, and.

38Suppose the domain Ω(t) has polygonal boundaries ΓI(t) and ΓO(t). Then we can cover the closure of by a regular triangulation of triangles, i.e., where each T is a closed triangle. The triangular mesh, which is a basic thing that Finite Elements requires, is generated by distmesh, which is a mesh generation tool implemented in MATLAB, and our algorithm is outlined in Algorithm S1. For the detailed implementations, such as: construct basis functions over the triangulation, assemble the stiffness matrix, etc, see references,.

39Numerical simulation is initialized by a small formed plaque. (see Figs. 4, 5, 6). Five combined levels of LDL and HDL ( and ) are tested for 300 days:

40: a small plaque doubles in size at 300 days;

41: a small plaque increases approximately 50% at 300 days;

42: a small plaque remains small at 300 days;

43: a small plaque decreases approximately 30% at 300 days;

44: a small plaque almost disappear at 300 days.

45Fig. 4 shows the growth of the plaque in case (a), Fig. 5 shows the shrinkage of the plaque in case (c), and Fig. 6 shows almost no plaque in case (e). In Fig. 7, the weight of the plaque, the summation of total cells, namely,, is plotted for these five scenarios of combined levels of LDL and HDL. Similarly to Fig. 7, we show in supporting information files how the populations of macrophages, SMCs, foam cells and T cells, as well as the concentration of ox-LDL, IFN-γ and IL-12, vary for different levels of LDL and HDL shown in Figs. S1, S2, S3, S4, S5, S6, S7.

46Fig. 8 shows a risk-map of plaque development. To create the risk-map, we divided the LDL axis by 121 equidistant points, i.e,, and divided the HDL axis by 21 equidistant points, i.e.,. For each pair, we computed the weight of the plaque, after 100 days on the domain shown in Fig. 3 (B), and formed the risk matrixwhere is the initial weight of the plaque. The vertical axis on the right of Fig. 8 shows the legend of the percentage of plaque growth or shrinkage. Accordingly, we divided the LDL-HDL plane into three regions: region I predicts high risk of plaque development, region III predicts no risk, and the intermediate region II predicts low risk.

47In order to support the robustness of the simulation results, we ran sensitivity analysis on parameters which appear in the differential equations and in the boundary conditions. The parameters chosen are those whose baseline was somewhat crudely estimated while at the same time they seem to play an important role in the development of the plaque. Specifically, we chose all the 15 production rate parameters from the third box of Table 4, all the 5 influx rate parameters from the third box of Table 5, and,. We list all these parameters with their range, baseline and unit in Table 6.

48Following the sensitivity analysis method described in, we performed Latin hypercube sampling and generated 100 samples to calculate the partial rank correlation coefficients (PRCC) and p-values with respect to the weight of the plaque after 300 days. The PRCCs are shown in Fig. 9, and all the p-values (not shown here) are less than 0.01. A positive PRCC (i.e., positive correlation) means that an increase in the parameter value will increase the weight of the plaque while a negative PRCC (i.e., negative correlation) means increase in the parameter will decrease the weight of the plaque. We note that is positively correlated, as it should be. Indeed, if is increased then MMP (Q) is increased (Eq. (12)) so that ECM (ρ) is decreased (Eq. (19)) and hence the plaque weight is increased (Eq. (16)). As another example, note that is negatively correlated. Indeed, if is increased then is decreased (Eq. (3)), and in the boundary condition will decrease, leading to smaller M, and then to smaller T and F. Similar explanation can be given to the other parameters.

49The most significant positively correlated parameters are and its influx rate. This is not surprising since LDL initiates the plaque formation. The most significant negatively correlated parameters are and its influx rate. Indeed, since HDL reduces the availability of free radicals, it plays an important negative role in plaque formation.

50Atherosclerosis is a disease in which a plaque builds up inside an artery. The process of plaque formation begins when, as a result of a lesion in the artery, cholesterols LDL and HDL enter the intima, and LDL becomes oxidized by free radicals. Upon sensing ox-LDL, endothelial cells secrete MCP-1 which attracts monocytes from the blood. As monocytes enter to the intima, they differentiate into macrophages that ingest the ox-LDL and become foam cells. Foam cells attract more macrophages, followed by T cells from the blood, and SMCs from the media. HDL reduces the available free radicals, as well as inflammation within the evolving plaque, thus HDL acts to block plaque growth.

51Public health guidelines in the U.S. specify that LDL level of 100–129 mg/dL is near ideal, 130–159 mg/dL is borderline high, and 160–189 mg/dL is very high, whereas concentration of HDL above 60 mg/dL is best, and below 40 mg/dL for men or below 50 mg/dL for women is poor. An important question is how to evaluate the risk of atherosclerosis for a pair of LDL and HDL taken together. This question is addressed in the present paper. We built a mathematical model of plaque development by a system of partial differential equations. The model includes two parameters:, the level of LDL in the blood, and, the the level of HDL in the blood.

52The model can simulate the evolution of a small plaque for any pair of values of. In Figs. 4, 5, 6, we simulated the plaque evolution over a period of 300 days. For example, one extreme case of mg/dL, mg/dL, the plaque doubled after 300 days; in another extreme case of mg/dL, mg/dL, the plaque disappeared after 300 days. We created a risk-map by taking sampling points of LDL and HDL values, and computing the weight of the plaque for each pair after 100 days. The map shown in Fig. 8, indicates the percentage of plaque growth or shrinkage for any such pair. We accordingly divided the (LDL,HDL) quadrant into three regions: high risk, low risk, and non risk.

53The need to consider the ratio of LDL/HDL in predicting coronary heart disease was suggested in a case study by. The American Heart Association considers the ratio to indicate high risk of heart disease, and the ratio to be risk free, where Tc denotes the total cholesterol, which is calculated by the formula Table 7 shows the National Cholesterol Education Program (NCEP) guidelines associated with plaque buildup. Accordingly, for the five points (a)–(e) in Results we have:

54for any value of Tr;

55if Tr is above normal;

56if Tr is not very high;

57if Tr is normal;

58if Tr is normal;.

59According to the NCEP guidelines, (a) and (b) should be in the high risk region; (c) in the low risk region; and (d), (e) in the no risk region, as indeed they are according placed in the risk map in Fig. 8.

60Some anti-cholesterol drugs, such as statins, lower LDL and at the same time also increase the HDL. It is important to know which drugs can best achieve the desired risk-free balance between LDL and HDL, that is, bring the individual's into the no risk (or low risk) region. By focusing not on just reducing or on just increasing, but on moving the combined to the no risk (or low risk) region in the shortest medically feasible path, we believe one could choose a more personalized medicine from those currently available, which will reduce the risk of atherosclerosis with the lowest amount of doze, thereby also possibly reducing potential negative side effects.

61To illustrate this approach, we note that for some drugs the ratio of decrease in LDL to increase in HDL is already known. For example, this ratio is 1/3 for the new experimental drug Evacetrapib. Some anti-cholesterol drugs only decrease LDL (e.g. Colestid) while others only increase HDL (e.g. Lofibra). We can the represent effect of such drugs by unit vectors in the plane: for example, Colestid ←, Lofibra ↑, and Evacetrapid. In Fig. 10. we consider three individuals, A, B, and C in the high risk region. In order to move them to the low risk region with the minimum amount of medication (side effects are ignored), the individual should choose the drug for which the line segment from the individual initial position to the low risk region is the shortest (We assume that the amount of drug is proportional to the length of the line segment). Thus the best drug for A is the one that primarily increases HDL. Similarly, C will do better with a drug that primarily decrease LDL, and B should use a drug with appropriate ratio of decreasing LDL to increasing HDL.

62Some work has been done on antioxidant therapy for reducing the risk of atherosclerosis, but so far it has had limited success in preventing cardiovascular diseases –. A review of studies in which antioxidant gene therapy has been successfully used is given in. Our model could account for antioxidative medication once we gain a good understanding of how such medication affects the source of free radicals, r 0, in Eq. (4).

63Some of the parameters in the differential equations in our model had to be rather crudely estimated, since no data were available, while others may slightly vary depending on the individual. As more data become available, parameter values may be further refined. Our model uses only the values of LDL and HDL as biomarkers. It will be interesting in the future to incorporate also triglycerides into the risk-map. Future work should also explore how other risk factors, such as high blood pressure, smoking and diabetes affect the risk-map.

64We did not include in this paper the circulation ox-LDL in the blood, which is elevated only in patients with advanced atherosclerosis,. Our model could be extended to include this additional biomarker, but at present there is not enough data on how the level of ox-LDL in the blood correlates to a specific advanced state of the disease.