Intrinsic Frequency Analysis and Fast Algorithms.

Intrinsic Frequency (IF) has recently been introduced as an ample signal processing method for analyzing carotid and aortic pulse pressure tracings. The IF method has also been introduced as an effective approach for the analysis of cardiovascular system dynamics. The physiological significance, convergence and accuracy of the IF algorithm has been established in prior works. In this paper, we show that the IF method could be derived by appropriate mathematical approximations from the Navier-Stokes and elasticity equations. We further introduce a fast algorithm for the IF method based on the mathematical analysis of this method. In particular, we demonstrate that the IF algorithm can be made faster, by a factor or more than 100 times, using a proper set of initial guesses based on the topology of the problem, fast analytical solution at each point iteration, and substituting the brute force algorithm with a pattern search method. Statistically, we observe that the algorithm presented in this article complies well with its brute-force counterpart. Furthermore, we will show that on a real dataset, the fast IF method can draw correlations between the extracted intrinsic frequency features and the infusion of certain drugs.

Introduction

Cardiovascular diseases (CVDs) relate to different conditions which can affect the performance of the heart and blood vessels. These diseases include: coronary artery disease, valvular heart disease, cardiomyopathy, heart rhythm disturbances and heart infections. As an example, congestive heart failure is a disease which happens when the heart does not work normally and cannot provide enough blood flow for body tissues. In this disease, the heart muscle does not stretch and contract in a normal way. Another example of CVDs is pulmonary hypertension which is a condition that there is a high blood pressure in the arteries that go from heart to the lungs. In some cases, the arteries in the lung become narrow or blocked. Therefore, blood flow harder through the arteries and this causes blood pressure to further increase in lungs making the blood flow even harder in this organ. All these changes will affect arterial pressure waveforms.

Cardiovascular diseases and stroke are major causes of death in the United States. The total cost related to CVDs and stroke was estimated to be more than $316 billion in 2011–2012[1,2]. The American Heart Association (AHA) estimates that in 2030 the direct costs will reach $818 billion[3]. It is mentioned that the indirect costs for 2030 would be $276 billion[3]. Hence, clinical measurements of cardiovascular health indices are of great importance. These methods and measurements are essential tools for monitoring cardiovascular health due to their relative availability. For example, Left Ventricular Ejection Fraction (LVEF) is a measure of left ventricular contractility[4]. Carotid-Femoral Pulse Wave Velocity (cfPWV) is a measure of aortic stiffness[5]. Velocity-encoded (VENC) phase contrast Magnetic Resonance Imaging (MRI) helps to detect the existence of flow vorticity that measures the morphological changes of the cardiac chamber wall[6].

However, current methods of measuring such indices are expensive, sometimes invasive, prone to measurement errors, and not necessarily easy to use. For example, 2D LVEF echocardiography is not accurate compared to more expensive and laborious gold standard cardiac MRI method[7-10]. As another example, obtaining accurate cfPWV measurements often requires certain medical devices and a well-trained staff within a clinical setting[11]. Consequently, continuous measurement of these indices is not practical. These limitations emphasize the need for new cardiovascular monitoring methods.

Intrinsic Frequency (IF) has been established as a new method of cardiovascular monitoring through a novel signal processing methodology[12]. The IF method needs only an uncalibrated pulse pressure[13] to extract pertinent information regarding the cardiovascular health of an individual[12]. The IF method has also been shown to be capable of non-invasively measuring LVEF by means of an iPhone camera[14]. We believe that methods like IF are of clinical and financial benefit in addressing cardiovascular monitoring.

In this paper, at first, we provide an overview of the IF method. Next, we present an approximate derivation of the IF model by combining Navier-Stokes equations and continuity with elasticity equations. This helps to build a solid mathematical foundation for the IF method and the analysis that follows. Later, we analyze the IF algorithm in the space of feasible solutions, and based on that, we introduce a new version of the IF algorithm which is faster than the current brute-force IF method[13] while maintaining the same accuracy. We then perform a case study on real pressure waveforms drawn from canine data using our new algorithm. We will see that the fast IF algorithm is capable of capturing the effects of different drug infusions on a canine subject.

Brief Overview of IF Method

A History of Analyzing Cardiovascular Pulse Waveform

Blood pressure was first measured by Hales in 1735[15]. In his measurements, he found that blood pressure is not constant in the arterial system. He related these variations to the elasticity of the arteries[15]. Currently, it is known that the shape of the arterial pulse wave is intimately related to the physiology and pathology of the whole arterial system[16]. There has been much research on analyzing the dynamics of blood pressure and flow in arterial systems[17-21]. Specifically, there are two main approaches to analyzing cardiovascular pulse wave data. One approach is based on a systematic mathematical framework for the cardiovascular system. The other is based on directly analyzing the pulse pressure waveform using signal processing methods.

An example of the systematic framework can be seen with the set of Windkessel models[22]. The formulation of a minimal lumped model of the arterial system was first presented by Westerhof et al.[22]. Based on a Windkessel model, the arterial system dynamics have been modeled through a combination of different elements such as resistance, compliance and impedance. In this simplified model of the arterial system, the blood flow dynamics is modeled by the interaction between the elements (assuming the blood flow acts as the current in the system). Because of the type of modeling, the wave transmission of the blood flow is neglected. As a result, the Windkessel models is not able to represent the entire dynamics of the blood flow in an arterial system accurately.

On the other hand, there are various methods for direct analysis of an arterial pulse waveform, in both time and frequency domains[20]. For example, the impedance method, which is based on Fourier transform, is a common method to analyze the pressure waveform in the frequency domain[17]. As an example, Milnor has shown that the pressure and flow waveforms can be a superposition of several harmonics using the Fourier method[23]. Another method to investigate the pressure wave in the time domain is the wave intensity analysis which is based on wavelet transform[24]. These methods do not necessarily convey a physical understanding of the cardiovascular system.

The IF algorithm presented in[13] is analyzing a pulse waveform through a direct time-frequency signal processing machinery setting, from a quantitative perspective. Although, in previous work[12], we tried to qualitatively express a systems approach to the IF formulation, the quantitative picture has not yet been clear. However, in this article, we show this connection from a quantitative perspective.

IF Formulation

In the IF method, the aortic pressure waveform at time , for a cardiac period T, can be represented aswith a continuity condition at T0 and periodicity at T. In this formulation, the indicator function is defined asHere, a1, b1, a2 and b2 are the envelopes of the IF model fit. ω1 and ω2 are the Intrinsic Frequencies (IFs) of the waveform. Further, is the mean pressure during the cardiac cycle. This type of formulation embeds the coupling and decoupling of heart and aorta.

The goal of the IF model (2.1) is to extract a fit, called Intrinsic Mode Function (IMF), that carries most of the energy (information) from an observed pressure waveform f(t) in one period. The latter is done by solving the following optimization problem[13]:subject toHere, is the L2-norm. The first linear condition in this optimization enforces the continuity of the extracted IMF at the dicrotic notch. The second one imposes the periodicity. The mathematical convergence and accuracy of the IF algorithm have been explained in a previous work[13]. In the next sections, we explore the foundation of the IF algorithm and propose a faster IF algorithm.

Approximate Derivation of the IF Model

As mentioned earlier, in a previous work[12], we tried to express a systems approach to the IF formulation qualitatively. However, in this article, we show this connection from a quantitative perspective. This section is devoted to this purpose.

Here we use a simplified model to address our approach. However, for a more general modeling, analysis and estimation of the blood flow and pressure estimation one could see[25,26].

In this paper, we assume that the Left Ventricle (LV), the aortic valve, aorta and the arterial system can be represented by a simplified model as shown in Fig. 1. Here, the LV and the aortic valve are assumed to be the boundary condition at the entrance of the aortic tube and the arterial system is the terminal boundary condition of the aortic tube. The boundary condition at the entrance of the aortic tube changes from an LV boundary condition to a closed valve boundary condition, at the dicrotic notch time T0 during a cardiac cycle [0, T]. We further assume that blood is a Newtonian incompressible fluid, the aorta is a straight and sufficiently long elastic tube with a constant circular cross section and there is no external force causing flow rotation. These assumptions are not all satisfied in a real cardiovascular system. However, they are useful in estimating the general behavior of blood in aorta.

Figure 1

Simplified cardiovascular system model schematic.

Combining the Navier-Stokes equations and continuity with the elasticity equation, we can drive a model for the flow Q(x, t) and the pressure P(x, t) along the length x of an aorta as follow

The step by step derivation of these equations is presented in the supplementary material. Parameters L, R, and C represent inductance, resistance, and compliance of the blood in aorta. Here, 0 ≤ x ≤ h, where h represents the aortic length. This model has also been discussed and simulated numerically in[27] with a complex set of boundary conditions. Here, our main concentration will be on the aortic tube oscillatory waveform solutions. Next, we will show that we can derive (2.1) from (3.1) and (3.2).

Since the input to the IF model (2.1) is a pressure waveform, we need to extract an equation for the pressure P(x, t) from Equations (3.1) and (3.2) by eliminating the flow. Combining Equations (3.1) and (3.2) results in

A similar expression could also be found for the flow field Q(x, t). Taking , with as the the mean pressure, we can write equation (3.3) asHere, we have used the dot notation to represent the time derivative. We can simplify the term in front of , in (3.4), by setting . The latter has a solution for some constant K. This reduces Equation (3.4) intoThe solution of Equation (3.5) can be expressed in terms of eigenfunctions. In other words, using the method of separation of the variables, one can express the solution of Equation (3.5) asforand some constants α, β, ζ and η. As a result, the solution of (3.3) can be expressed as

The variables ω can be expressed based on the boundary conditions of the aortic tube. We need to emphasize that for a period of the cardiac cycle [0, T), the boundary conditions change before and after the dicrotic notch T0. Hence, for , Equation (3.9) can be written asHere, the superscripts indicated with “1” belong to the form of the solution before the closure of the aortic valve, and the superscripts indicated with “2” belong to the form of the solution after the closure of the aortic valve. Constants K1, , , , and are found from the boundary and initial conditions at systole. Similarly, constants K2, , , , and are found from the boundary and initial conditions at diastole.

Equation (3.10) is explicitly showing the coupling and decoupling of heart and aorta before and after the dicrotic notch. As the boundary conditions change during a cardiac cycle, the frequencies of oscillation also change from to . Generally, Equation (3.9) can represent pressure waveform for a Newtonian incompressible fluid in a straight and sufficiently long elastic tube with constant circular cross section.

If the pressure is recorded at a specific point x0 on aorta, the terms containing the spacial variable x would be fixed. In other words, Equation (3.10) would reduce toforand

Now, considering that the cardiac cycle length would be around 1.5 sec, at most, and taking into account that R is smaller than L[27], one can use the approximation . Hence, Equation (3.11) would become

Further, if most of the information, or energy, is carried out by the first terms in the series of the solution, we can further write the approximated solution (3.12) as

Now, by relabelingwe can approximate the IF model (2.1). The continuity and periodicity conditions (2.3) can also be approximated if we hold the assumption that most of the energy is carried out by the first terms in the series of the solution (3.12).

In short, in this section, we have presented an approximate quantitative justification on the origins of the IF method. In the next section, we move on with the analysis of the optimization problem (2.2) subject to (2.3).

Analysis of The IF Algorithm

Practically, one must solve the discrete version of (2.2). We assume that the pressure waveform f(t) is sampled uniformly. Also, we can simplify (2.2) by the fact that any sinusoid can be assumed to start from time t = 0 with a compensation coming from a phase shift. In other words, any sinusoid can be expressed as A cos ωt + B sin ωt, irrespective of whether the initial time is t = 0 or t = T0. Hence, the discrete format of (2.2) can be expressed asfor f = (f1, …, f)′ as the uniform sampling of the original cycle. Here, by takingfor and , we have the discrete form of as

In this article, (.)′ denotes the transpose operator and the vector . Also,

The constraints, in (4.1), can be written as

If we can solve for two, out of four, unknowns in (4.5), we would make (4.1) an unconstrained optimization. However, it is important to check whether the matrix in (4.5) is of full rank or not. In fact, the rows of this matrix are linearly independent except when

This will lead into two cases:

Degenerate Case in which Equation (4.6) holds,

General Case in which, it does not.

General Case (cos ω1T0 cos ω2(T − T0) ≠ 1)

One can solve the constraints in (4.1) for a1 and a2 to obtainEquations (4.7) and (4.8) would then simplify (4.3) intowhere Q(ω1, ω1, b1, b2; t) = b1v1(ω1, ω2; t) + b2v2(ω1, ω2; t) forandUsing Equations (4.9–4.11), and dropping the dependencies in notation, simplifies (4.1) into

This simplification has helped to eliminate the constraints in the optimization problem (4.1).

The minimization problem (4.12) is non-convex and non-linear in its parameters. So, in order to be able to solve the problem, we can use the fact that the minimum of a function can first be found over some variables and then over the remaining ones[28]. In other words, the optimization problem in (4.12) can be written as

We call the inner optimization in (4.13) as P(ω1, ω2). Solving for P(ω1, ω2) is a classical least squares problem. The solution existence and uniqueness of this optimization is mentioned in our previous work[13]. To find the exact solution we simplify the objective function as

Substituting for Q = b1v1 + b2v2, we convert (4.14) into

Since, in this part of the optimization, the values of ω1 and ω2 are fixed, we can find the optimal values of b1, b2, and by setting the partial derivatives of (4.15) equal to zero. In other words, we set , , and . Doing this, we find the optimal solution for b1, b2, and , byHere, we have fulfilled the optimization part by solving a linear system. This could potentially accelerate the IF algorithm. Finally, we only have to solve a minimization onwhich is

We note that a property of the function P(ω1, ω2) is its differentiability, away from its singularities. In fact, by definition, the function is directionally differentiable with respect to all its variables. Hence, using the results in[29,30], we can deduce thatis directionally differentiable with respect to ω1 and ω2. This property can be exploited if one tries to solve (4.18) using a gradient based optimization method[31].

Degenerate Case (cos ω1T0 cos ω2(T − T0) = 1)

The solution of (4.6) can be expressed as nodes of a lattice in ω1ω2 plane. To be more specific, we havewhereand

If (ω1, ω2) ∈ Γ1, from (4.5) we have a1 = −a2. On the other hand, if (ω1, ω2) ∈ Γ2, from (4.5) we have a1 = a2. In both of these cases, we can express (4.3) aswhere   =  , for i = 1, 2. If (ω1, ω2) ∈ Γ1,

Similarly, if (ω1, ω2) ∈ Γ2, we have

In both of the cases, we haveandHere, 01 and 02 are zero vectors in and , respectively. It is clear, from (4.26) and (4.27), that . Using (4.23), and a similar approach we employed in (4.15) and (4.16), we find the optimal solution for a1, b1, b2, and , byfor i = 1, 2. Hence, similar to (4.17), for (ω1, ω2) ∈ Γ1 or (ω1, ω2) ∈ Γ2, we only have to solve a minimization on

Note that, from a machine learning perspective, the nodes specified in (4.20) do not have important information physiologically as they could be inferred from the systolic and diastolic parts of a waveform alone. In other words, even if these points present a global minima, they are not informative as we already know the systolic and diastolic inverses, and respectively, as possible inputs to any machine learning algorithm. Hence, these points could possibly be ignored in a search for an optimum point of (4.1).

Fast IF Algorithms

In this section, we present a fast IF algorithm which is based on the results presented in the previous section and the topology of the solution space for P(ω1, ω2). In order to keep the fluency of this section, we mention the original IF algorithm (see Algorithm 1) as presented in[13].

Algorithm 1 has three major steps. In the first step, the (ω1,ω2) domainis made discrete, namely . The second step is a minimization to find P(ω1, ω2), see (4.18). The final step is a brute-force search on to find the minimum of P(ω1, ω2).

All three steps can be optimized to make the IF algorithm faster. Regarding the domain of optimization , defined in (5.1), we know from our previous work in[14] that the average IF solution, for a physiological pulse waveform recording, is confined to a smaller domain expressed as

This will make the first step search area more well-defined and optimized. In the previous section, we have been able to find some analytic solutions (see (4.16)) for the inner optimization part of problem (4.13). This will help us to substitute an analytic solution instead of an iterative[32] or QR decomposition[33] solution for (5.3). Finally, the brute-force part can be substituted with an appropriate direct search algorithm[34], e.g. pattern search algorithm[35]. It can even be substituted with an appropriate gradient based algorithm[28,31], e.g. gradient descent, as we know the differentiability of P(ω1, ω2).

Intrinsic Frequency.

Before moving on, we show the topology of the P(ω1, ω2) function and also its minima locations in ω1 and ω2 space. These will provide useful insights on where to set the initialization point(s) of a possible fast IF algorithm. The data description is provided in the next section. In Figs 2 and 3, we have presented two different dog aortic pressure cycles with the IMF extracted by the means of the brute-force IF Algorithm 1. Figures 2 and 3, top right, show the heat-map plots of . The complex nature of P(ω1, ω2) can be seen in these figures. We purposefully plotted P in the dimensionless coordinates and to show the behavior of this function with respect to the lattice node locations defined in (4.20–4.22). To have a better view and understanding of the P(ω1,ω2) topology, a contour of is shown in those figures. The general topology of , for all aortic or carotid pulse waveforms, is similar to the ones presented in Figs 2 and 3. However, the location of the minimizer is not similar.

Figure 2

Up-Left: A dog aortic pressure cycle (in blue), its dicrotic notch (black circle with dotted center), and the IMF (in red). Up-Right: heat-map plot of for the cycle in left with the location of the solution marked with red dot. Down-Left: Contour plot of . The location of the minimizer of P is shown by a red dot. Down-Right: Contour plot of and the location of the minimizer of P tracked by the pattern search. The beginning of the pattern search is marked with a triangle and its end with a circle. The true optimum point is marked with a red dot. The upper pattern search set (in black) has converged towards the correct optimum. The lower pattern search set (in magenta) has converged to a local minima near the node.

Figure 3

Up-Left: A dog aortic pressure cycle (in blue), its dicrotic notch (black circle with dotted center), and the IMF (in red). Up-Right: heat-map plot of for the cycle in left with the location of the solution marked with red dot. Down-Left: Contour plot of . The location of the minimizer of P is shown by a red dot. Down-Right: Contour plot of and the location of the minimizer of P tracked by the pattern search. The beginning of the pattern search is marked with a triangle and its end with a circle. The true optimum point is marked with a red dot. The lower pattern search set (in black) has converged towards the correct optimum. The upper pattern search set (in magenta) has converged to a local minima near the node.

Our investigations show that the locations of the minimizers of all P functions construct two different areas in the dimensionless coordinates and . We call these areas as the upper lobe and lower lobe. The upper lobe is an area, in , confined above the line . The lower lobe is an area, in , confined below the line . This is also the case for human subject data[14]. This type of topology suggests two critical initial guess areas for any non-brute-force algorithm solving (4.1): one set of points in the upper lobe, the other in the lower. In the remaining part of this section, we introduce a fast IF algorithm based on the pattern search method[34].

Pattern Search IF

The pattern search algorithm (or sometimes called the compass search algorithm) is explained in detail in[34]. For completeness, we have summarized the pattern search algorithm in Algorithm 2. The convergence analysis of this method is expressed in[34].

Pattern Search[34].

The fast IF algorithm, without considering the nodes (4.20), is expressed in Algorithm 3. As mentioned before, what makes Algorithm 3 fast is embedded in three different objects:

The initial guess set up in the initialization part of the algorithm.

The fast analytic solution at each point iteration defined by (4.16) and (4.17).

The pattern search part which is a substitute for the brute force algorithm.

Figure 2, bottom right, shows the results of Algorithm 3. In this figure, when using Algorithm 3, we have used two initial guesses and . As depicted on the figure, the initial guess located in the upper lobe has converged towards the true minimizer in . On a PC having 8 threads, Intel® Core™ i7-4700MQ CPU @ 2.40 GHz × 8, running a Matlab implementation of the brute-force Algorithm 1 in parallel takes roughly 85 seconds. On the other hand, achieving the same minimizer, using a sequential version of the fast Algorithm 3, takes approximately 0.5 seconds. Comparing the brute-force and pattern search algorithms, the absolute error in estimating ω1 is ~0.02 (an absolute error of ~0.002 in ) and the relative error is 0.17%. The absolute error in estimating ω2 is ~0.08 (an absolute error of ~0.008 in ) and the relative error is 0.74%.

The same test was done for another aortic cycle presented in Fig. 3. We used the same initial guesses as before. This time, on the same PC, using the same implementations, the brute-force Algorithm 1 took roughly 80 seconds and the fast Algorithm 3 took approximately 0.5 seconds. These two examples show a speed up of almost 160 times. In the next section we present more about the statistical accuracy of Algorithm 3 and its physiological capabilities. Comparing the brute-force and pattern search algorithms, the absolute error in estimating ω1 is ~0.02 (an absolute error of ~0.001 in ) and the relative error is 0.11%. The absolute error in estimating ω2 is ~0.05 (an absolute error of ~0.006 in ) and the relative error is 0.83%.

Fast IF.

Real Data Example

Data Description

The real dog data used in this manuscript is well described in[36]. Briefly, six normal adult beagle dogs had undergone the data collection experiment. One dog was involved in a sterile surgical procedure for implanting chronic recording transducers. All the dogs had undergone general anesthesia with propofol and maintained it with inhaled isoflorane. For the aortic pressure waveform measurements, a micro-manometer-tipped catheter was inserted into a femoral artery and guided into the descending thoracic aorta. The transducer outputs were transfered to a personal computer through an A/D conversion system. The measurements were collected over a period of 50– 170 minutes. Some pharmacological influences and total intra-vascular volume changes were imposed on the dogs: dobutamine, esmolol, verapamil, phenylephrine, nitroprusside, saline and progressive hemorrhage. For a detailed description, see[36].

Since, at the time of the the data retrieval, the data was downloaded with different sampling rates, we re-sampled all six dog data at 500 Hz. We used a modified version of the automatic cycle selection introduced in[37] to pick cycles. Dicrotic notch locations were then found from the picked cycles[38]. We totally extracted 59384 acceptable aortic cycles form those six dogs.

Statistical Accuracy

To check the statistical accuracy of the fast IF algorithm versus the brute-force IF algorithm, we compared the results of these two algorithms on the extracted 59384 dog aortic cycles. The brute-force IF algorithm (Algorithm 1) was run over the sample set with a mesh size of 0.02π. Algorithm 3 was run on the same sample set of 59384 aortic cycles with Δω = 0.001, and Δω0 = 0.1, comprising a mesh size of . The brute-force algorithm has a larger mesh size due to heavy computational cost of this algorithm. The maximum average difference between the IFs found by these two algorithms was found to be less than 0.0475. This difference is smaller than both mesh sizes used for the brute-force and fast IF algorithms. This shows that, on average, the fast IF algorithm (Algorithm 3) reaches the same minima as the brute-force algorithm (Algorithm 1).

Physiological Observations

To evaluate the new fast IF algorithm (Algorithm 3), we applied the algorithm on the measured aortic pressure signal from one dog experiencing various pharmacological interventions, see Fig. 4. During the experiment, the dog was under the following pharmacological influences: infusion of dobutamine (5–20 μg/kg/min), phenylephrine (2–8 μg/kg/min) and nitroglycerin (4 μg/kg/min) during different time intervals.

Figure 4

Drug effects on ω1 and ω2. First Panel: The measured aortic pressure waveform recorded in time. Second Panel: Dosage of dobutamine (in red), phenylephrine (in purple), and nitroglycerin (in green) during the aortic pressure measurement. Third Panel: Changes of ω1 in units of bit per minute (bpm) over the measurement time. Each drug effect is projected with its corresponding color, for the fast IF algorithm. No drug areas are in blue, for the fast IF algorithm. The brute-force algorithm results are shown in dashed black line. Fourth Panel: Changes of ω2 in units of bpm over the measurement time. Each drug effect is projected with its corresponding color, for the fast IF algorithm. No drug areas are in blue, for the fast IF algorithm. The brute-force algorithm results are shown in dashed black line.

The second panel, in Fig. 4, shows the dosage and duration of each drug in the experiment. In the first phase of the experiment dobutamine has been injected at a low dosage followed by a fluctuation in the dosage of injection.

The third and fourth panels, in Fig. 4, show the trends of ω1 and ω2 for both the fast IF algorithm (Algorithm 3 in blue, red, purple and green dots) and the brute-force algorithm (Algorithm 1 in dashed black line). Both algorithms follow the same trends showing the same precision. In other words, if the force-brute IF is used as a control, the fast IF is exactly duplicating the trend.

The effect of dobutamine on the cardiovascular system is to increase the strength and force of the heartbeat. Consequently, it forces more blood to circulate throughout the body. In previous works[12,14], we hypothesized that ω1 would be a representative of heart functionality. We also hypothesized that ω1 and ω2 would try to keep a balance during changes. These hypotheses can be seen during the injection of dobutamine in this figure.

Next, phenylephrine has been injected at a low dosage and the dosage is then increased over time. Phenylephrine is a decongestant which affects the cardiovascular system by shrinking blood vessels. ω2 shows an almost monotone decrease during the infusion of phenylephrine. This is again in qualitative accord with what we presented in[12].

Lastly, nitroglycerin has been injected at a constant dosage. Nitroglycerin helps to dilate the blood vessels. This dilation can be captured with ω2, as can be seen from the figure. Generally, based on this figure, IFs are able to capture changes in the dynamics of the system under the effects of different drugs.

Conclusion

In this paper, we provided a mathematical foundation for the IF model[13]. We showed how to derive an estimation of the IF model (2.1) by considering basic physics principles. More precisely, we showed that the IF model can be estimated from Navier-Stokes and elasticity equations.

We further analysed the IF model (4.1). This helped to introduce a fast algorithm for the IF method (Algorithm 3). What made this algorithm fast was embedded in the proper set up of the initial guesses based on the topology of the problem, fast analytic solution at each point iteration, and substituting the brute force algorithm with a pattern search method. These changes would convert an iterative and brute-force method (Algorithm 1) into an algebraic and iterative method (Algorithm 3). The presented fast algorithm, in this article, has a speed up of more than 100 times compared to the brute-force algorithm provided in[13]. From a statistical perspective, we have also shown that the algorithm presented in this article complies well with the brute-force implementations of this method.

We also showed, on a real dataset, that the fast IF Algorithm 3 can depict correlations between its outputs and infusion of certain drugs. This part of our paper can be subject to further physiological and clinical investigations in a future work.

Data availability

All data generated or analysed during this study are included in the Supplementary Information files. The datasets generated during and analysed during the current study are also available from the corresponding author on reasonable request.

Research ethics

All experiments and procedures were reviewed and approved by the MSU All-University Committee on Animal Use and Care[36].