A Particle Filter Based Framework for the Prognosis of Atherosclerosis via Lumped Cardiovascular Modeling

Atherosclerosis refers to the plaque deposition in the arteries that can eventually lead to any of the three cardiovascular diseases, namely, heart attack, stroke, or peripheral vascular disease, depending upon the site of the blockage in the human arterial network. This work attempts to prognose this pathological condition via lumped cardiovascular modeling while utilizing the radial artery blood pressure measurements. To achieve this, the cardiovascular system has been modeled as a third order non-linear system with explicit emphasis on the systemic circulation. The parameters of the model are estimated using non-linear least squares estimation technique by minimizing the error between the measured and the estimated arterial pressure waveforms. The arterial pressure is found to be sensitive to three of the model parameters, namely, arterial compliance, systemic vascular resistance, and the peak cardiac muscle elastance. Based on the analysis, a growth model of systolic blood pressure is developed as a function of the arterial blockage. A particle filter based mathematical framework is then utilized to predict the time it would take to reach the stage of critical arterial blockage that may cause heart attacks.


INTRODUCTION
Cardiovascular diseases (CVDs) are the leading cause of the deaths across the globe as per World Health Organization.An estimated 17.7 million people died from CVDs in 2015 that accounts for 31 % of all global deaths.It can be noted that most of the CVDs can be prevented by addressing behavioral risk factors such as unhealthy diet and obesity, tobacco usage, harmful alcohol usage, and physical inactivity.But once the CVD has been diagnosed, it needs immediate treatment, otherwise it may prove to be fatal.brovascular disease (stroke), peripheral vascular disease, heart failure, rheumatic heart disease, congenital heart disease and cardiomyopathies.In this work, we shall focus on the category of cardiovascular diseases that arise because of the arterial blockage (usually referred to as atherosclerosis).Depending upon the site of the blockage in the human arterial network, atherosclerosis can lead to three types of CVDs.If the arteries supplying blood to the brain are blocked, the pathological condition is referred to as stroke.If the arteries supplying blood to the heart are blocked, the condition is referred to as coronary artery disease.If the arteries supplying blood to the lower extremities of the body are blocked, the condition is referred to as peripheral vascular disease.In this work, we shall be focused on the cause, that is, atherosclerosis.Later in this paper, we shall be using the phrase "cardiac health status (CHS)" as an indicator of the arterial blockage.
Near the beginning of the 20th century, German physiologist Otto Frank formulated the two-element Windkessel model that was viewed as the load against which the heart had to pump the blood; resistor and capacitor were the two elements there.The resistor represented the resistance of the arteries and the arterioles, together known as the total peripheral resistance.The capacitor represented the capacitance of the large arteries, better known as the total arterial compliance.Later on, modifications of this model were proposed, namely, three-element and four-element Windkessel models (Westerhof, Lankhaar, & Westerhof, 2009).Though wave transmission and wave travel could not be studied via these 'lumped' models; their simplicity and the ease of implementation made them a preferred choice over distributed models to understand the pressure-flow relations at the entrance of the arterial system, which is of clinical relevance.
The instantaneous left ventricular pressure-volume ratio was shown to be independent of mechanical loading conditions; however it changed with changes in contractile state of the ventricle (Suga & Sagawa, 1974).This ratio is referred to as cardiac muscle elastance and it has long been used to model the dynamics of the atrial and/or ventricular chambers of the heart.The ventricle described using time-varying elastance function is coupled to the Windkessel arterial load; this network is filled from a constant venous pressure reservoir (Yu, Boston, Simaan, & Antaki, 1998).In the last six decades, many models at different levels of complexity have been reported in the literature to describe the cardiovascular dynamics (Williams et al., 2014), (Jain, Maka, & Patra, 2018).While the higher order models can be potentially more accurate and physically meaningful, the major pitfall of the higher order models is the large number of parameters associated with them.It becomes quite difficult to estimate all the parameters with the available measurements, which are often obtained under controlled situations.The proposed work utilizes a third order subject-dependent non-linear model to describe the cardiovascular dynamics.
One of the earlier works has developed a mathematical model of plaque formation via partial differential equations with the aim to determine the risk of plaque formation for combined levels of low density lipoprotein (LDL) and high density lipoprotein (HDL) (Hao & Friedman, 2014).Their work creates a "risk map" for plaque formation in the LDL-HDL coordinate plane.In recent works, the prognostic values of coronary and systemic atherosclerosis in the conditions of coronary artery disease and myocardial infarction respectively have been investigated (Assante et al., 2017), (Calais et al., 2018).Most of the prognostic works in this direction have been carried out while utilizing statistical models or distributed models.To the best of our knowledge, the proposed work is a first attempt to utilize a lumped cardiovascular model for the prognosis of atherosclerosis.
In order to overcome the aforementioned risk factors associated with the human heart, a prognostic framework can be very useful in detecting the systolic blood pressure (SBP) growth from an early stage of life and in providing a sufficiently accurate estimate of the cardiac health status (CHS) of a person suffering from atherosclerosis.Based on the literature it can be observed that prognostics study has been carried out by researchers on numerous applications such as structural damage prognosis (Farrar & Lieven, 2007), crack growth analysis (Sankararaman, Ling, Shantz, & Mahadevan, 2011), Li-ion batteries (Guha & Patra, 2018), human neuromusculoskeletal system (Mussleman, Gates, & Djurdjanovic, 2016), etc. Regarding approaches to prognostics study, various methods are available in the literature such as Autoregressive Integrated Moving Average (ARIMA), Extended Kalman Filter (EKF) and Particle Filter (PF) (Saha, Goebel, & Christophersen, 2009).As ARIMA is an exclusively data-driven approach, it does not consider any physics of the system during computation, thereby ending up with wide uncertainty margins.Hence it is not a good choice for long term predictions.In case of EKF, although it is quite robust against nonstationarity, it is not capable of accommodating un-modeled effects and thus may diverge quickly.Whereas a PF not only provides a mean estimate of the time-to-failure but also gen-erates a probability distribution over time that best characterizes the uncertainties within the system model and measurements.In the work by (Guha & Patra, 2018), internal resistance growth model of a lithium-ion battery, based on electrochemical impedance spectroscopy (EIS) test data, is developed to mark the aging of the battery.This is analogous to the growth in blood pressure of a human with the gradual occlusion of the arteries.In the battery model, parameters are initialized using a curve fitting approach; they are then updated with the new measurements while utilizing a particle filter (PF) framework.Tuned parameters are then utilized to predict the remaining useful life (RUL) of the battery.In our work, we have adopted a similar approach to predict the time that would elapse before atherosclerosis reaches a critical stage where heart attacks may occur.
The remaining paper is organized as follows.Section 2 covers the methodology adopted for the cardiovascular assessment.Subsection 2.1 explains the electrical circuit based model of the cardiovascular system.It also includes the mathematical framework of non-linear least squares (NLLS) estimation routine adopted to estimate the model parameters.Subsection 2.2 shows how we generated synthetic systolic blood pressure (SBP) data from the model parameters.Subsection 2.3 explains how SBP gets updated via particle filtering.Subsection 2.4 computes the time required for the blockage to reach a critical level.Section 3 highlights the important results and is divided into two subsections.Subsection 3.1 presents the results from the parameter estimation part.Subsection 3.2 shows results from the particle filtering approach that has been utilized to predict the time it would take to reach the stage of critical arterial blockage that may cause heart attacks.Section 4 concludes the paper.

Modeling and Estimation of Cardiovascular System
Blood circulation in the body can be viewed as an electrical system in which the heart acts as a voltage source and the rest of the body tissues form the systemic load.The volume of blood in the left ventricular chamber of the heart, the blood pressures in the systemic arteries, and veins form the three state variables of the model (Jain, Patra, & Maka, 2018).The equivalent electrical circuit for blood circulation is presented in Figure 1.The model is based upon pressure-voltage analogy.Resistors here represent resistance to the blood flow.Capacitors indicate the compliance of the vessels.Blood flow is analogous to current.Pressure in the vessels is analogous to the capacitor voltage.The complete list of the physiological variables and the parameters is presented in the Nomenclature section.
The ordinary differential equations representing the blood circulation are derived using circuit theory concepts.The state vector for the model is defined as follows: Figure 1.The proposed cardiovascular circulation model T Consider the node marked as v h in the circuit diagram.The flow entering this node would be equal to the flow leaving this node.That is, the flow through the mitral valve resistance R mv would be equal to the rate of change of the ventricular volume v h plus the flow through the aortic valve resistance R av .Mathematically, This can be re-written as: This is further explained as follows.The ventricular dynamics are truly described by the time-varying compliance which is just the reciprocal of the elastance function E h (t) described in the paper.However, from circuit theory viewpoint, we presented the heart's ventricular chamber as a voltage source representing a pump, as it pumps the blood from low (venous) pressure side to high (arterial) pressure side.The rest of the equations, derived using the nodal analysis, are listed as follows. (2) where, The cardiac muscle elastance is quantitatively expressed as follows, (Williams et al., 2014).
The following two expressions are used for the valvular resistances. (5) where, the open valve resistance, R op is 0.001 units and closed valve resistance, R cl is 20 units.
The main assumptions considered while modeling are: • The right side of the heart as well as the pulmonary circulation have been lumped along with the veins; only the systemic side of the circulation has been considered explicitly.
• The large vessels, namely, the systemic arteries and veins have been approximated by their compliances whereas the small vessels, namely, peripheral portion of the circulation have been approximated via systemic vascular resistances.
Furthermore, the problem of minimizing the error between the model output, namely, the arterial pressure and the actual blood pressure data leads to the tuning of the model parameters, thereby, making the model subject-dependent.The following cost function is minimized to obtain the optimized set of parameters: where, p m a,sys,l and p c a,sys,l are respectively the measured and the computed values of systolic arterial blood pressure in l th cardiac cycle; p m a,dia,l and p c a,dia,l are respectively the measured and the computed values of diastolic arterial blood pressure in l th cardiac cycle, and M denotes the number of cardiac cycles.It can be noted that classical sensitivity analysis is used to find out the relatively more sensitive parameters (Batzel, Bachar, & Kappel, 2012).Only these parameters are tuned and the remaining ones are kept at their initialized values.The set of relatively more sensitive parameters for the proposed model found out to be:

Systolic Blood Pressure (SBP) Data Generation
Let us recall Eq. ( 8) that contains most sensitive parameters w.r.t the arterial pressure p a .This vector θ opt is utilized to generate the mean values of systolic blood pressure (SBP) at different instants (months).The need to generate this data is explained as follows.We have the arterial blood pressure data collected for a finite number of seconds.From this data, we can compute the present mean value of the systolic component of the blood pressure.We need to predict how this mean value would vary with time considering the present lifestyle of an individual.If an individual has very healthy lifestyle, it is assumed that the mean value of SBP would not rise rapidly as the time progresses.However, if the individual has a sedentary lifestyle, the mean value of SBP is assumed to rise rapidly.This much of information is utilized to decide the rate of change of the parameter vector θ opt .The vector changes at a higher rate for an individual maintaining a sedentary lifestyle.
Before we proceed further, it is worth mentioning that an arterial blockage is often related to high blood pressure (WebMD, 2016).It can be noted that the systolic blood pressure of 180 mmHg can be considered as the threshold value (American Heart Association, 2014).Therefore it can be said that the arterial blockage becomes critical when the pressure reaches its threshold value.Now, we shall be discussing the data generation in more detail.
Five young healthy subjects each having 20 years of age, namely, H1, H2, H3, H4, and H5 are considered and their arterial blood pressure along with the demographical information is recorded.The subjects are deliberately chosen to be of the same age so that a comparative prognostics study could be done while considering their present lifestyle.It would be shown via a growth model of SBP that a subject maintaining healthy lifestyle in terms of diet and exercise would reach the threshold SBP value in greater time span as compared to those having a sedentary lifestyle.Now, let us consider subject H1, having gender, height and weight as male, 180 cm and 75 kg, respectively.The optimized set of parameters are obtained as per the procedure explained earlier.Let us again recall the set of sensitive parameters (with respect to the arterial pressure) denoted by θ opt . Here, This parameter vector (excluding E min , which is comparatively less sensitive) is varied gradually in successive steps to obtain the arterial pressure, p a at each of those steps.The mean value of the systolic blood pressure is recorded at these steps.Let us assume that H1 maintains a healthy life style and therefore the parameter variation r i (see Figure 3) is assumed to be as low as 1 % at every step.These steps are mapped into the number of months.The blockage increases to its threshold value in greater time span as compared to other subjects.It can be noted that the parameter C a is increased, and the parameters R p and E max get reduced by the same proportion.It is explained as follows.As the plaque builds up, compliance of the vessel gets reduced; resistance to the blood flow increases and therefore cardiac muscle contracts with greater force to push the same amount of blood to the systemic tissues.We further performed this numerical experiment for four more healthy subjects, namely, H2, H3, H4, and H5.They are assumed to have the parameter variation at the rates of 1.5, 2.5, 3.0, and 3.5 respectively; all values are in percentages.The higher the rate of variation, the poorer the life style in context to health.This means that the person with the highest rate has a sedentary life style; does not keep healthy diet and avoids exercise, and thereby catches up with the threshold blockage in the shortest time span.
While utilizing the parameter vector θ opt , the regression model capable of providing a good fit to the systolic blood pressure values at different instants is obtained.The regression equation is given by: where k is a time instant (in months), and p sys,k is the systolic blood pressure at the k th instant.The fitted curve for the representative subject is shown in Figure 2. We shall now proceed to develop the mathematical framework for the prognostics approach.

SBP Model Updating by Particle Filter Algorithm
Dynamic state estimation and prediction do not depend only on an accurate model but also rely on model parameter tuning to track the variation in health degradation.In order to identify the dynamic health degradation, the particle filtering (PF) approach is used to estimate the current health param-eter.The basic idea of PF is based on recursive Bayesian filtering using Monte Carlo simulations.In the PF framework (Arulampalam, Maskell, Gordon, & Clapp, 2002), the probability density function (pdf) is approximated by a set of particles constituting sampled values from an unknown state space along with a set of associated weights, thereby indicating discrete probability masses.Particles are generated and recursively updated based on a probabilistic model in addition to the available measurements.
The coefficients p 1 -p 4 of the regression model (see Eq. ( 9)) form the state vector s of the particle filter.The mathematical framework for the filter is presented as follows.The state equations are as follows: where, j = 1, ..., 4, and q j,k ∼ N (0, σ 2 pj ).The measurement function z (= SBP ) is as follows: Here, q j,k is the process noise, and vn k is the measurement noise.Both the process and measurement noises are assumed to follow a normal distribution with zero mean and variances σ 2 pj (=1e-08) and σ 2 SBP (=0.16) respectively.The weights are such that they sum up to unity: where, w i k denotes one particle, and N p represents the total number of particles.The particles approximate the distribution of s k by: where, s i k denotes a single estimate of the state s k .With the advancement of k, the posterior pdf becomes the prior pdf of the new k with the dynamic model f (.).Here, P (s k |z k ) denotes the posterior pdf and P (s k |z k−1 ) denotes the prior pdf.
Once the prior distribution of s k is available, the posterior distribution of s k can be obtained by: where, P (z k |s i k ) denotes the likelihood of s i k .The likelihood is calculated by: where, σ represents the standard deviation of the measurement noise vn k .In order to satisfy the condition that the weights sum up to unity as given in Eq. ( 12), normalization is performed as follows.
The w i k in Eq. ( 17) denotes the normalized weight.

Prediction of Cardiac Health Status (CHS)
While utilizing the regression model of Eq. ( 9), the h step ahead prediction (He, Williard, Osterman, & Pecht, 2011) of each trajectory at instant k is obtained as: The estimated posterior pdf of the prediction can be obtained for each trajectory along with the associated weights as follows: The mean value of the h step ahead prediction at the k th instant can be obtained as: Based on the assumption of the failure threshold limit as α times the systolic blood pressure at the healthy state of the subject, the CHS estimation C i k of the i th trajectory of the k th instant can be computed by: where p sys,th indicates threshold value of SBP, and α is defined as the ratio of the threshold value of SBP to its initial stage value p sys,healthy .
Therefore, C i k is the value of the CHS corresponding to the i th trajectory when it crosses the threshold p sys,th .After obtaining all these C i k values for i = 1 to N p , the CHS pdf is computed as follows: Basically it is the posterior estimate of the probability density function provided the measurement is available up to the time instant k.For each particle i, a trajectory is obtained which intersects the threshold line at some point of time.Similarly, for N p particles, N p such trajectories are obtained which intersect the threshold line at different time instants.Finally, based on the mean and variance of all those intersecting points on the threshold line, the CHS pdf has been obtained in this work.This concept is numerically illustrated in the Results section.The mean value of the CHS prediction at the k th instant is given by: It may be noted that we shall be using the abbreviation SBP for the systolic blood pressure p sys , for convenience.
To summarize, an overall block diagram illustrating the proposed approach for the prognosis of atherosclerosis (CHS prediction) is presented in Figure 3. Firstly, the optimized parameter vector θ opt is obtained for a given subject using nonlinear least squares (NLLS) algorithm.Secondly, the parameter vector is varied at a degradation rate r k .The corresponding mean values of SBP are recorded at different instants.Next, the regression model for the SBP growth is obtained.The coefficients of this model are then updated via a particle filter.Finally, the cardiac health status (CHS) is predicted.In the block diagram, θ 0 denotes set of initialized pa-Figure 3. The proposed approach for the prognosis of atherosclerosis (CHS prediction) rameters, θ opt denotes optimized parameters, r k is the subject dependent degradation rate and it equals the rate at which the parameters would change every instant k, the generated parameters would be referred as θ k , SBP k is the corresponding value of systolic blood pressure, n is Gaussian noise such that n ∼ N (0, 0.16), and N p indicates number of particles which is numerically considered to be 100.

Modeling and Parameter Estimation
Let us consider the aforementioned representative subject, H1.The model parameters are initialized using empirical formulas that require subject's age, gender, height, weight along with the blood pressure recording (Williams et al., 2014).Further, the classical sensitivity analysis revealed that four of the model parameters are relatively more sensitive, namely, arterial compliance, C a , peripheral resistance, R p , and peak elastance values, E max and E min .Therefore, only these parameters are optimized using non-linear least squares estimation technique (as described in Section 2); the remaining ones are kept at their initialized values.The arterial blood pressure waveforms for the subject H1, before and after optimization, are shown in Figure 4.The corresponding parameter values

Prognosis of Atherosclerosis via Particle Filter Framework
In order to analyze the proposed method, the error corresponding to the SBP threshold (p sys,th ) along with its standard deviation have been considered as two indicators in terms of the number of months to assess the prediction accuracy.
A lower value of the standard deviation indicates a prediction with a higher confidence interval.The prediction result corresponding to the 285th month is shown in Figure 5a where only the data corresponding to the first 285 months have been considered to update the model.Similarly, the prediction results corresponding to the 420th and the 555th months have been shown in Figures 5b and 5c respectively.The presumed threshold value for the growth in SBP with respect to a healthy subject, as shown in Figure 5, is chosen to be 180 mmHg.
From Table 2, it can be observed that with the advancement of the prediction time towards the threshold value, the prediction accuracy improves gradually.It is observed that the subject H1 who leads a healthy lifestyle (as mentioned earlier) would reach the threshold value of the critical arterial blockage at the end of 996 months, that is, at the 83 years of age.On the other hand, the subject H5 that leads a sedentary lifestyle (as mentioned earlier) would probably reach the threshold value of the arterial blockage as early as at the age of 564 months, that is, 47 years.In the

CONCLUSION
In the paper, the human cardiovascular system is modeled as a third order lumped parameter model with explicit emphasis on the systemic circulation.The model parameters are estimated using non-linear least squares estimation technique by minimizing the error between the measured and the estimated arterial pressure waveforms.A growth model for systolic blood pressure is developed that relies on three of the The dynamic model (in context of PF) g(.) The

Figure 2 .
Figure 2. The fitted SBP curve for the subject H1 Figure 4. (a) Systemic arterial pressure waveforms-measurement and the estimated, before optimization; (b) the waveforms after optimization

Figure 5 .
Figure 5. CHS prediction results based on SBP growth of subject H1 at 285th, 420th and 555th months respectively

Table 1 .
(Williams et al., 2014)erse problem converges to one local minimum(Jain, Patra, & Maka, 2018).It may be further noted that the healthy blood pressure data has been taken from open resources made available by North Carolina State University(Williams et al., 2014).

Table ,
observation function (in context of PF) Diastolic value of p a (computed) in l th cycle p m a,dia,l Diastolic value of p a (measured) in l th cycle p c a,sys,l Systolic value of p a (computed) in l th cycle p m a,sys,l Systolic value of p a (measured) in l th cycle p sysThe h step ahead prediction of i th trajectory of p sys at k th instant p j j th coefficient of SBP regression model P (.)