A mathematical model to assess the effects of COVID-19 on the cardiocirculatory system | Scientific Reports – Nature.com

The modified lumped-parameter model consists of a system of ordinary differential equations (ODEs) that needs to be numerically solved to allow the computation of different model outputs of clinical interest. We calibrated the model to fit some clinical data of patients hospitalized for severe COVID-19-related pneumonia in the Internal Medicine ward of L. Sacco Hospital in Milan, Italy, between March and April 2020. We analysed the statistical reliability of the model outputs for each successful calibration by means of uncertainty intervals and, finally, we performed a statistical analysis on clinical data or model outputs by means of hypothesis tests to highlight the impairments of the cardiocirculatory system associated with COVID-19 pneumonia.

We identified four groups of quantities, taken from the dataset or obtained as an output of the calibrated model:

The clinical data used for the model calibration, obtained from clinical measurements and referring to physical quantities (PQ1), as, for example, the maximal left atrial volume (LAVmax) and the systolic systemic pressure (SAPmax);

The inputs of the model (heart rate HR and body surface area BSA) and of the calibration procedure (right ventricular fractional area change RVFAC and tricuspid annular plane systolic excursion TAPSE), provided by other clinical measurements;

The parameters of the model (e.g. resistances and compliances) determined through a calibration procedure, from now on referred to as calibrated parameters;

The outputs of the numerical simulation of the model (e.g. flow rate and mean pressure), from now on referred to as model outputs. Some of them (MO1) referred to physical quantities (PQ1) that were also measured (clinical data), for example, LAVmax and SAPmax. Other model outputs (MO2) referred to physical quantities (PQ2) that were not measured but quantified only by means of the computational model. Examples of the latter are the mean left atrial pressure (LAPmean) and indexed right ventricular end diastolic volume (RVI-EDV). The complete list of PQ1 and PQ2 is reported in Supplementary Table 1.

We remark that the indexed value of volumes of a patient can be computed dividing the volumes by the BSA of that patient (Supplementary Table 2). In what follows, an I- that precedes a subscript of a volume means that the volume is indexed (for example, LVI-EDV is the indexed left ventricular end diastolic volume).

For the sake of clarity, we reported in Fig.1 the diagram flowchart of the followed procedure that is described in detail in what follows.

Diagram flowchart of the procedure used in this study. Top: calibration; mid: statistical analysis of measured physical quantities; bottom: statistical analysis of computed physical quantities. Calibration: the mathematical model required as inputs HR and BSA of a specific patient. The model computed MO1 using an initial setting of parameters (that could need to be calibrated, so they are highlighted in red). If MO1 were close enough to the clinical data the model was considered calibrated (the parameters are highlighted in green); if not, the calibration method was iteratively applied to the parameters using RVFAC and TAPSE as inputs. If the parameters were not modified the calibration failed; if not, MO1 were recomputed by using the new setting of parameters and the previous steps were repeated. Statistical analysis 1: we performed hypothesis tests on clinical data (test I). Statistical analysis 2: HR and BSA were used as inputs of the calibrated model for every patient with a successful calibration, the model computed the MO2 and we checked the statistical reliability of MO2. We collected the reliable MO2 from every patient and we performed hypothesis tests on the reliable MO2 of all the patients (test II).

The dataset consists of (58) patients, who all required oxygen supplementation but none of them was on mechanical ventilation. Of such patients, only (29) were calibrated according to point (iii) above (see Calibration subsection below) ((56pm 18) years). Such patients did not present symptoms or signs of heart failure or substantial structural cardiac disease; (10) out of (29) were older than (64) years; (6) patients had arterial hypertension, (1) had diabetes and (4) showed the association of hypertension and diabetes (Supplementary Table 2).

The echocardiography of each patient was performed early after the admission to the hospital. Examinations were performed at bedside using a Philips CX-50 portable device by expert operators. Measures were defined according to the latest European and American Echocardiography Society guidelines18,19.

Each patient provided consent to use his/her data for observational studies. The institutional board has approved the study with protocol number 16088/2020.

The cardiovascular system was studied by means of a lumped-parameter (0D) mathematical model that splits the system into compartments (e.g. right atrium, systemic arteries/veins) and, for each of them, the time evolution of model outputs (pressures, flow rates and cardiac volumes) is modelled by a system of ODEs20,21. The lumped-parameter model is described through an electrical circuit analogy: the current represents the blood flow through vessels and valves; the electric potential the blood pressure; the electric resistance plays the role of the resistance to blood flow; the capacitance represents the vessel compliance; the inductance the blood inertia; the increase in elastance the cardiac contractility.

There are different possible choices and number of compartments, depending on the purpose of the study, for the construction of a lumped-parameter model (e.g. Refs.16,17,22,23). We considered the computational model introduced in Ref.17, wherein the four heart chambers, the systemic and pulmonary circulations, with their arterial and venous compartments were included, and we substituted the 3D left ventricle with a 0D component (as in Ref.16) and we added two new compartments accounting for systemic and pulmonary capillaries. The pulmonary capillary circulation was also split in two compartments accounting for oxygenated and non-oxygenated capillaries (Fig.2).

Lumped-parameter cardiocirculatory model. The unknown pressures and flow rates are in red and blue, respectively, whereas the model parameters are in black. Notice in the green boxes the new compartments with respect to16 featuring this work.

The system of ODEs associated with the lumped-parameter model is formed by the equations representing continuity of flow rates at nodes and of pressures in the compartments, and its numerical solution allows to compute several model outputs as functions of time: the left and right atrial and ventricular volumes (({V}_{{text{LA}}}), ({V}_{{text{LV}}}), ({V}_{{text{RA}}}) and ({V}_{{text{RV}}})), the systemic and pulmonary arterial, capillary and venous pressures (({p}_{{text{AR}}}^{{text{SYS}}}), ({p}_{{text{C}}}^{{text{SYS}}}), ({p}_{{text{VEN}}}^{{text{SYS}}}), ({p}_{{text{AR}}}^{{text{PUL}}}), ({p}_{{text{C}}}^{{text{PUL}}}) and ({p}_{{text{VEN}}}^{{text{PUL}}})), the systemic and pulmonary arterial and venous blood fluxes (({Q}_{{text{AR}}}^{{text{SYS}}}), ({Q}_{{text{VEN}}}^{{text{SYS}}}), ({Q}_{{text{AR}}}^{{text{PUL}}}) and ({Q}_{{text{VEN}}}^{{text{PUL}}})).

Starting from these functions, it is possible to compute the pressures of the four cardiac chambers (({p}_{{text{LA}}}), ({p}_{{text{LV}}}), ({p}_{{text{RA}}}) and ({p}_{{text{RV}}})), the blood fluxes through the valves (({Q}_{{text{MV}}}), ({Q}_{{text{AV}}}), ({Q}_{{text{TV}}}) and ({Q}_{{text{PV}}})), through the systemic capillaries (({Q}_{{text{C}}}^{{text{SYS}}})) and through oxygenated and non-oxygenated pulmonary capillaries (({Q}_{{text{C}}}^{{text{PUL}}}) and ({Q}_{{text{SH}}})), and all the model outputs referring to PQ1 and PQ2 (Supplementary Table S1).

We considered reference values of the parameters (such as resistances and compliances) such that all the model outputs were in the reference healthy ranges of the corresponding physical quantities taken from the literature7,18,19,24 for an ideal individual with HR equal to (80) bpm (beats per minute) and BSA equal to (1.79) m2 (Supplementary Table S3). We did not consider model outputs computed starting from the flow rates, because they are not uniquely defined depending on the tract of the compartment where they are measured, from ({p}_{{text{C}}}^{{text{SYS}}}), due to the heterogeneity of the pressures of systemic capillaries among tissues, and from ({p}_{{text{VEN}}}^{{text{SYS}}}), even if we recovered the value of central venous pressure, that coincides with the right atrial pressure24.

We reported the system of ODEs associated with the lumped-parameter model in Supplementary Equations S1. The lumped-parameter model was numerically discretized by means of Dormand-Prince method25 (adaptive stepsize RungeKutta) which was implemented in Python using the Jax library26.

The lumped-parameter model was characterized by parameters representing the functional properties of the compartments (e.g. resistances). To properly select such values for a specific compartment and patient, a calibration procedure was needed27,28.

We chose a priori the cardiac timings and the resistance of oxygenated pulmonary capillaries (({{text{R}}}_{{text{C}}}^{{text{PUL}}})) equal to the associated reference values. In particular, we fixed ({{text{R}}}_{{text{C}}}^{{text{PUL}}}) to avoid modelling micro-thrombosis because of its possible increase. For the remaining parameters, the calibration of the model relied on the method we presented in Ref.27, that is aimed to reduce the sum of squared relative errors between the model outputs MO1 and clinical data, modifying the parameters of the model in suitable bounded intervals ({{text{I}}}_{{text{i}}}), for (i=1,dots ,{N}_{{text{p}}}), where ({N}_{{text{p}}}) is the number of parameters, independent of the patient, built starting from the reference values of parameters mentioned before (Supplementary Table S3). Specifically, we chose to calibrate those parameters among the latter according to a sensitivity analysis estimating the absolute correlation coefficients between parameters and model outputs (Supplementary Table S4). We calibrated only the parameters featuring at least one absolute correlation coefficient greater than (0.1) that was associated to provided clinical data. To reproduce the blunted hypoxic pulmonary vasoconstriction condition, the resistance of non-oxygenated pulmonary capillaries (({{text{R}}}_{{text{SH}}})) could decrease in such a way that the shunt fraction could reach values up to (70%) in the worst-case scenario. The list of amendable parameters varies between different patients according to the different clinical data provided.

The calibration was based on clinical measurements of COVID-19 patients that were provided by L. Sacco Hospital in Milan and referred to HR and BSA, which were used as inputs for the lumped-parameter model, RVFAC and TAPSE, which determined the bounded interval ({{text{I}}}_{overline{{text{i}}} }) used during the calibration, with (overline{i }) the index referring to the right ventricular active elastance, and the clinical data, given by a subset of the pressures and volumes involved in the cardiac circulation (Supplementary Table S2).

To provide further mathematical details, we indicate with (mathbf{p}) a configuration of parameters of the cardiocirculatory model. The calibration method aimed to find the configuration of parameters ({overline{mathbf{p}} }^{{text{j}}}) which minimized the loss function for the specific patient (j), that reads:

$$L^{{text{j}}} left( {mathbf{p}} right) = sumlimits_{{{text{l}} = 1}}^{{{text{N}}^{{text{j}}} }} {left( {frac{{q_{{{text{m}}_{{text{j}}} left( {text{l}} right)}}^{{text{j}}} left( {mathbf{p}} right) - d_{{text{l}}}^{{text{j}}} }}{{d_{{text{l}}}^{{text{j}}} }}} right)^{2} } ,$$

(1)

where ({N}^{{text{j}}}) is the number of available echographic clinical data for patient (j), ({d}_{{text{l}}}^{{text{j}}}) is the value of the l-th clinical data of patient (j) (Supplementary Table S2) and ({q}_{{{text{m}}}_{{text{j}}}left({text{l}}right)}^{{text{j}}}) is the value of the model output related to the l-th clinical data of patient (j). The index (m) of ({q}_{{text{m}}}^{{text{j}}}) lies in ({1,dots ,{N}_{{text{q}}}}) where ({N}_{{text{q}}}) is the number of both MO1 and MO2. We considered the model calibrated for a specific patient if the loss function was below ({10}^{-3}). Notice that, for some patients, the calibration procedure could fail, if, for example, it reaches the minimum of the loss function that is above the required threshold.

Moreover, to improve the robustness of the calibration procedure, we repeated, for every patient, the calibration three times, with different initial configurations of parameters, and we considered the calibrated setting of parameters that returned the lowest loss function. As anticipated above, only (29) out of (58) patients were successfully calibrated. We noticed that by performing (4) times the calibration procedure the number of calibrated patients was still equal to (29), precisely as after (3) calibrations.

The loss function (1) was minimized by the Quasi-Newton method L-BFGS-B29 implemented in Scipy by computing its gradient by means of automatic differentiation (reverse mode gradient) included in the library Jax26.

For every patient (j) calibrated with a loss function below ({10}^{-3}), a configuration of parameters ({overline{mathbf{p}} }^{{text{j}}}) was at disposal. The loss function was computed using the clinical data provided by L. Sacco Hospital, which were related to measurement errors (Supplementary Table S1), that also affected the uncertainty of the model outputs ({mathbf{q}}^{{text{j}}}). We needed to determine, for every patient, if the related model outputs were reliable or not, so we proceeded along two steps:

Build a sample of candidate model outputs ({mathbf{q}}^{{text{j}},{text{k}}}) for (k=1,...,n) ((n) was (100));

Determine, by employing a simple statistical analysis, whether the mean of the model outputs was reliable.

Regarding step 1, for every provided clinical data ({d}_{{text{l}}}^{{text{j}}}) of patient (j), we built an interval ({{text{M}}}_{{text{l}}}^{{text{j}}}) centred in the value of the clinical data with width equal to two times the measurement error (Supplementary Table S1). Then, we built the samples ({mathbf{q}}^{{text{j}},{text{k}}}) by following the subsequent procedure:

Choose a relative width (w) ((w) was (12.5%));

Build an interval centred at ({overline{p} }_{{text{i}}}^{{text{j}}}) and with width (2w{overline{p} }_{{text{i}}}^{{text{j}}}) for every (i=1,dots ,{N}_{{text{p}}}). If this interval is not included in the parameter interval ({{text{I}}}_{{text{i}}}) used for the calibration, then cut off its overflowing extremities.

Perturb every parameter of the calibrated patient sampling from a uniform distribution in the corresponding interval built at point b) thus obtaining ({p}_{{text{i}}}^{{text{j}}});

Run a simulation of the cardiocirculatory model with parameters ({mathbf{p}}^{{text{j}}});

Check if the model output ({q}_{{{text{m}}}_{{text{j}}}left({text{l}}right)}^{{text{j}}}) generated at point d) lie in the intervals ({{text{M}}}_{{text{l}}}^{{text{j}}}). If they do, save the new configuration of acceptable model outputs ({mathbf{q}}^{{text{j}}}), otherwise reject it;

Repeat from point c) until (n) iterations are performed;

Check if the acceptance ratio (ratio between the number of saved configurations and the number of iterations) is within ([0.1, 0.15]). If it does, repeat from point c) to e) until (n) configurations are accepted because at this step the sample size of candidate model outputs is small (with (n=100), the size is between (10) and (15)), otherwise increase or decrease (w) to retrieve the condition on the acceptance ratio, discard the previous configurations and repeat from point b).

Once the above procedure was concluded, we proceeded with step 2 by using the (n) samples of acceptable model outputs ({mathbf{q}}^{{text{j}},{text{k}}}) for (k = 1, dots , n) generated at the previous step, for every specific patient (j). If the standard deviation of the sample of a model output of patient (j) was lower than 5% of its mean, we considered the mean reliable and we used it for the hypothesis tests. In this way, for every model output we built a sample of accepted values (depending on the patient), where sample size depended on the considered model output.

Prediction intervals could have been used for this analysis, but, if the sample was not normally distributed, a link function would be needed to retrieve normality30. We checked, for every patient (j) and for every model output, if the sample of that model output was normally distributed by means of a chi-squared test. It turned out that the sample is not normally distributed for all patients. Thus, since we wanted to use the same statistical approach for every patient, we resorted to this heuristic approach based on standard deviation instead of prediction intervals.

If the sample mean, calculated over all patients, of a clinical data or MO2 (referring to physical quantities PQ1 and PQ2, respectively) fell inside the healthy range of the corresponding physical quantity7,18,19,24, we did not consider the physical quantity altered in association with COVID-19 infection, otherwise we performed hypothesis tests to check whether the mean was significantly (p-value below (0.01)) increased or decreased with respect to the healthy range to investigate the impairments of the cardiovascular system in association with COVID-19 infection. If the sample mean, calculated over all patients, was less than the lower bound of the healthy range, the null hypothesis was that the mean was greater or equal than the lower bound of the healthy range, whereas the alternative hypothesis was that the mean was smaller than the lower bound of the healthy range. If we accepted the null hypothesis, then the corresponding physical quantity was considered not altered in association with the infection of COVID-19; otherwise, we considered the physical quantity altered in association with COVID-19. If, instead, the sample mean was greater than the upper bound of the healthy range, we proceeded similarly.

For each clinical datum, we computed the mean and the standard deviation of its sample without resorting to the mathematical model. The sample sizes were large enough to use one-tailed z-tests (assuming the variance equal to the unbiased sample variance) comparing their means to the nearest bound of the healthy range (test I).

For every MO2 we computed the mean and the standard deviation of its sample. We performed a chi-squared test and not every sample was normally distributed, so we opted for one-tailed z-tests (assuming the variance equal to the unbiased sample variance) only if the sample had more than (24) elements comparing their means to the nearest bound of the healthy range (test II).

Notice that for group PQ1 the statistical analysis was carried out directly using the clinical data and not the MO1 values. Accordingly, the clinical data were used in a twofold way:

To statistically compare PQ1 clinical measures with healthy ranges independently of the application of the proposed lumped-parameter model (test I);

To calibrate the lumped-parameter model for the patients at hand thus allowing to obtain MO2 that are statistically compared with healthy ranges (test II).

See the original post here:

A mathematical model to assess the effects of COVID-19 on the cardiocirculatory system | Scientific Reports - Nature.com