Changes in Resting and Exercise Hemodynamics Early After Heart Transplantation: A Simulation Perspective

Introduction: During heart transplantation (HTx), cardiac denervation is inevitable, thus typically resulting in chronic resting tachycardia and chronotropic incompetence with possible consequences in patient quality of life and clinical outcomes. To this date, knowledge of hemodynamic changes early after HTx is still incomplete. This study aims at providing a model-based description of the complex hemodynamic changes at rest and during exercise in HTx recipients (HTxRs). Materials and Methods: A numerical model of early HTxRs is developed that integrates intrinsic and autonomic heart rate (HR) control into a lumped-parameter cardiovascular system model. Intrinsic HR control is realized by a single-cell sinoatrial (SA) node model. Autonomic HR control is governed by aortic baroreflex and pulmonary stretch reflex and modulates SA node activity through neurotransmitter release. The model is tuned based on published clinical data of 15 studies. Simulations of rest and exercise are performed to study hemodynamic changes associated with HTxRs. Results: Simulations of HTxRs at rest predict a substantially increased HR [93.8 vs. 69.5 beats/min (bpm)] due to vagal denervation while maintaining normal cardiac output (CO) (5.2 vs. 5.6 L/min) through a reduction in stroke volume (SV) (55.4 vs. 82 mL). Simulations of exercise predict markedly reduced peak CO (13 vs. 19.8 L/min) primarily resulting from diminished peak HRs (133.9 vs. 169 bpm) and reduced ventricular contractility. Yet, the model results show that HTxRs can maintain normal CO for low- to medium-intensity exercise by increased SV augmentation through the Frank–Starling mechanism. Conclusion: Relevant hemodynamic changes occur after HTx. Simulations suggest that (1) increased resting HRs solely result from the absence of vagal tone; (2) chronotropic incompetence is the main limiting factor of exercise capacity whereby peripheral factors play a secondary role; and (3) despite the diminished exercise capacity, HTxRs can compensate chronotropic incompetence by a preload-mediated increase in SV augmentation and thus maintain normal CO in low- to medium-intensity exercise.


INTRODUCTION
Heart transplantation (HTx) is the last resort for an increasing number of persons suffering from end-stage heart failure. During HTx, sympathetic and vagal denervation of the heart is inevitable, leading to postoperative chronotropic incompetence and thus to reduced quality of life (Awad et al., 2016). Especially in the first year after HTx, patients suffer from chronic tachycardia, with resting heart rates (HRs) elevated to greater than 90 beats/min (bpm) and significantly reduced HR variability (HRV) (Awad et al., 2016;Kobashigawa and Olymbios, 2017). Furthermore, on top of raised HR recovery times, HTx recipients (HTxRs) show delayed and impeded exercise response (Awad et al., 2016;Kobashigawa and Olymbios, 2017), reaching peak HRs as low as only 133 bpm (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Wilson et al., 1991;Marzo et al., 1992;Rudas et al., 1993;Kao et al., 1994;Doering et al., 1996;Geny et al., 1996;Notarius et al., 1998;Hayman et al., 2010;Peled et al., 2017;Nygaard et al., 2019;Nytrøen et al., 2019). According to literature, elevated resting HRs are most certainly a result of the absence of vagal tone, rendering the heart to rely on intrinsic control only, whereas delayed and impeded exercise response is most likely due to lack of vagal withdrawal and missing sympathetic drive (Awad et al., 2016). The exercise response of HTxR is thought to depend solely upon circulating catecholamines and a strong reliance on the Frank-Starling mechanism to compensate for the impaired autonomic response (Awad et al., 2016). Ultimately, elevated resting HRs to greater than 90 bpm and increased postexercise HR recovery times were shown to be strongly correlated with raised mortality (Awad et al., 2016;Kobashigawa and Olymbios, 2017).
With globally increasing prevalence, cardiovascular diseases are still the leading cause of mortality, projected to be responsible for more than 25 million global deaths by 2030 (Okwuosa et al., 2016;Timmis et al., 2020). This situation combined with increased donor organ availability leads to a steady increase of currently 5,500 annual HTx worldwide (Khush et al., 2019). Consequently, there is a major need for treatment modalities that help to improve the quality of life and clinical outcomes early after HTx (Awad et al., 2016). However, to this date, still, there is a substantial lack of knowledge on the mechanisms associated with cardiac denervation and spontaneous reinnervation in HTxRs (Awad et al., 2016), thus restraining the development of new therapies to provide relief for those affected. Therefore, tools and frameworks that facilitate the investigation of physiological changes in HTxRs to find potential treatment modalities are of constant need and great significance.
In recent years, numerical simulations have gained increasing popularity in cardiovascular research, especially to gain a better understanding of the human hemodynamics and its changes associated with various diseases (Jung et al., 2006;Banerjee et al., 2008;Qian et al., 2010;Politi et al., 2016). Moreover, mathematical modeling was successfully proven to be an effective technique for the design and evaluation of potential treatment modalities for heart failure patients, such as mechanical circulatory assist devices (Moscato et al., 2010(Moscato et al., , 2013Gross et al., 2020).
Although there are various numerical attempts to study the cardiac autonomic neural regulation (Levy and Zieske, 1969;Katona et al., 1970;Ursino and Magosso, 2003;Kember et al., 2014Kember et al., , 2017, to the best of our knowledge, up until now, there is no literature concerned with the simulation of hemodynamic changes and alterations of autonomic cardiac control associated with HTxRs. Therefore, based on published clinical data, we developed an integrated numerical model that is capable of predicting resting and exercise hemodynamics of HTxRs with good accuracy. The model does not attempt to represent a single patient's behavior but rather cohort hemodynamics.

MATERIALS AND METHODS
We present an integrated numerical model that was modified to reproduce the changes in hemodynamic parameters observed in early HTxRs, both at rest and during exercise. The proposed model integrates three closely linked components: (1) a closedloop lumped parameter model of the human cardiovascular system; (2) the intrinsic HR control represented by a Hodgkin-Huxley-type single-cell sinoatrial (SA) node model; and (3) the autonomic cardiac control mediated by arterial baroreflex and the pulmonary stretch reflex. A schematic overview of the model structure can be found in Figure 1. Finally, the model parameters were adjusted based on clinical data from 15 publications (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Wilson et al., 1991;Marzo et al., 1992;Rudas et al., 1993;Kao et al., 1994;Doering et al., 1996;Geny et al., 1996;Notarius et al., 1998;Hayman et al., 2010;Peled et al., 2017;Nygaard et al., 2019;Nytrøen et al., 2019) that were identified in a careful literature review. The model equations were implemented in MathWorks R SIMULINK R 9.1 (2019a) and numerically integrated using the ode15s solver with a maximum step size of 0.01 s, and an absolute error tolerance of 10 −3 . All parameter values used to simulate HTxRs and an age-and gender-matched healthy control group may be found in Supplementary Tables S2-S8.

Cardiovascular System
To simulate the human cardiovascular system, the numerical closed-loop lumped-parameter model published by Moscato et al. (2010Moscato et al. ( , 2013 was modified to reproduce the hemodynamics of early HTxRs and demographically matched healthy cohorts reported in the literature (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Wilson et al., 1991;Marzo et al., 1992;Rudas et al., 1993;Kao et al., 1994;Doering et al., 1996;Geny et al., 1996;Notarius et al., 1998;Hayman et al., 2010;Peled et al., 2017;Nygaard et al., 2019;Nytrøen et al., 2019). The hemodynamic model comprises atria and ventricles represented by time-varying non-linear elastance models, heart valves realized by diode resistance series, and compartments for systemic and pulmonary circulation with their venous and arterial portions. The temporal effects of respiration on intrathoracic and abdominal pressures were integrated into the model according to Ursino and Magosso (2003). FIGURE 1 | Schematic overview of the integrated numerical model, comprising (1) the cardiovascular system and its principal compartments; (2) the single-cell SA node model; and (3) the autonomic control governed by arterial baroreflex and pulmonary stretch reflex integrating receptor information on instantaneous aortic pressure (P ao ) and instantaneous lung volume (V L ). Based on autonomic activity, the SA node depolarization rate is modulated through changes in acetylcholine (ACh) and isoprenaline (Iso) concentrations. Autonomic control has also influence on the left and right ventricular contractility (E lv , E rv ), systemic vascular resistance (R as ), and systemic venous unstressed volume (V usv ). Muscle mechanoreflex and metaboreflex modulate Ras to account for vasodilation in response to exercise. Exercise intensity determines the level of central command, respiration frequency, and muscle mechanoreflex and metaboreflex activity. Ultimately, the central command has a direct influence on autonomic control loops and the adrenal medulla releasing catecholamines (Iso circ ).
The equations of the cardiovascular system may be found in the Supplementary Material (S1-S26), and the interested reader can find detailed information in Ursino and Magosso (2003) and Moscato et al. (2010Moscato et al. ( , 2013.

Intrinsic Heart Rate Control
The intrinsic HR control was realized by a Hodgkin-Huxleytype single-cell human SA node model (Pohl et al., 2016). To cover the physiological range of HR, gating variables of the hyperpolarization-activated current, L-type Ca 2+ current, acetylcholine (ACh)-activated K + current, Na + /K + pump current, slowly activating delayed rectifier K + current, and sarcoendoplasmic reticulum Ca 2+ -ATPase (SERCA) activity were adjusted according to Fabbri et al. (2017). To simulate the effects of time-varying neurotransmitter concentrations, gating variable equations were further adjusted according to Zhang et al. (2002) and Li et al. (2014). Thus, the cell model incorporates the influence of the parasympathetic neurotransmitter ACh, and the sympathetic neurotransmitters epinephrine and norepinephrine using the β-sympathomimetic drug isoprenaline (Iso) as a surrogate. The modified model produced a depolarization rate of 74 bpm, which corresponds to the values reported in Fabbri et al. (2017) and is also in good accordance with the experimental findings of Verkerk et al. (2007), who report an experimentally determined spontaneous beating rate of 73 ± 3 bpm for isolated human SA node cells.
The cell model equations are given in the Supplementary Material (S57-S175). Detailed information on the development of the SA cell model may be found in Zhang et al. (2002); Li et al. (2014) and Fabbri et al. (2017).
The instantaneous HR is derived from the time differences between the peaks of consecutive SA node action potentials and fed to the hemodynamic model.

Autonomic Control
The autonomic control of total peripheral resistance, venous unstressed volume, left and right ventricular elastance, and HR was modeled according to Ursino and Magosso (2003). However, the proposed sigmoidal function to calculate the HR in the original model (Ursino and Magosso, 2003) was replaced by the previously described single-cell SA node model in which the depolarization rate is determined by time-varying neurotransmitter concentrations. Furthermore, inspired by the work of Magosso and Ursino (2002), additive terms were introduced in all effectors accounting for the influence of central command, which contributes significantly to the exercise response. Finally, the autonomic control model was extended by two linear equations, responsible for the modulation of aortic and arterial compliance based on exercise intensity. To implement the HR modulation by neurotransmitter released from sympathetic and parasympathetic cardiac nerve terminals, sympathetic outputs were linearly correlated with changes in Iso, and vagal outputs with changes in ACh concentrations. Finally, the calculated concentrations were added to the respective baseline values to find the total neurotransmitter concentrations.
The fundamental principle of exercise response was implemented according to the considerations of Moscato et al. (2013) and therein cited literature (Rowell, 1993), including the interaction of central command and muscle mechanoreflex and metaboreflex. Additionally, the release of circulating catecholamines by the adrenal medulla in response to physical activity was introduced into the model as a first-order dynamic system. The direct effects of the central command were modeled as previously described by Magosso and Ursino (2002).
Finally, a gradual exercise intensity-dependent decrease of the respiratory period was implemented, starting at a basal resting period of 5 s, eventually reaching 1.2 s at peak exercise. To simulate the muscle mechanoreflex and metaboreflex, peripheral resistance was decreased by 50% from an initial value of 0.97 mmHg·s/mL. The operating point of arterial baroreflex was increased by 90% from a basal value of 91 mmHg.
All equations governing the autonomic control can be found in the Supplementary Material (S16-S46).

Parameterization to Simulate HTxRs
Heart transplantation leads to severe changes in the autonomic feedback loops of HR and ventricular elastance control. This can be mainly attributed to efferent cardiac denervation impeding the vagal and sympathetic neural outflow to the allograft. On top of that, HTxRs also exhibit peripheral cardiovascular abnormalities, including increased total peripheral resistance and reduced vasodilatory capacity.
In Ursino and Magosso's model of cardiac autonomic control (Ursino and Magosso, 2003), the magnitude of efferent outflow to the heart is determined by its cardiac effector gains. Therefore, to simulate cardiac denervation following HTx, modifications were applied to the gains of the heart period and ventricular elastance effectors. The vagal and sympathetic heart period effector gains were reduced by 100 and 95%, respectively. Furthermore, the sympathetic gain of ventricular elastance control was reduced by 95%.
To account for the peripheral changes in HTxRs, selected parameters of the hemodynamic model described by Moscato et al. (2010Moscato et al. ( , 2013 were modified as follows. The total peripheral resistance was increased by 15%, whereas the peripheral compliance and the vasodilatory effect of muscle mechanoreflex and metaboreflex were decreased by 20 and 10%, respectively. Possible changes in afterload resulting from HTx medication were implicitly incorporated into the model by tuning its parameters for pooled published data of HTxRs receiving typical post-HTx medication regimen (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Wilson et al., 1991;Marzo et al., 1992;Rudas et al., 1993;Kao et al., 1994;Doering et al., 1996;Geny et al., 1996;Notarius et al., 1998;Hayman et al., 2010;Peled et al., 2017;Nygaard et al., 2019;Nytrøen et al., 2019). A breakdown of the prescribed pharmacological agents by study can be found in Supplementary Table S1.
In all studies investigating exercise response (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Marzo et al., 1992;Kao et al., 1994;Notarius et al., 1998), modified square wave endurance exercise tests were performed in which the exercise intensity was increased until the maximum tolerated power (MTP) was reached. In one study (McLaughlin et al., 1978), the exercise test was performed in a supine position using an ergometer. In all other publications, the exercise tests were performed in upright positions, using either ergometers or treadmills.
Patients of all studies received a typical post-HTx medication regimen. In one publication (McLaughlin et al., 1978), no information on medication was given; therefore, we assumed conventional post-HTx medication as in the other studies. The breakdown of prescribed pharmacological agents by study can be found in Supplementary Table S1.
In case that the body mass index was not given, it was calculated from height and weight. In one study (Labovitz et al., 1989), data were given as median and range; therefore, mean was substituted by the median, and variance was calculated dividing the given range by four.
The results for pooled hemodynamic parameters can be found in Table 1.

Validation Population
For model validation, the simulation results for rest and peak exercise hemodynamics of HTxRs and healthy individuals were compared to published clinical data of an age-and gendermatched cohort other than that used for modeling but with comparable demographics (Pflugfelder et al., 1988;Hosenpud et al., 1989;Tamburino et al., 1989;Niset et al., 1991;Braith et al., 1992;Scott et al., 1995;Nytrøen et al., 2011;Cattermole et al., 2017) (Table 2). Because of a lack of data availability, the values used for validation of healthy resting hemodynamics were taken from a study reporting normative ranges of an agematched healthy cohort (Cattermole et al., 2017). The literature was selected from an online search on Google Scholar, using the previously described set of keywords and exclusion criteria.
The results for pooled hemodynamic parameters can be found in Table 2.
conditions. The protocols were designed with two goals in mind.
(1) The exercise protocol should be in accordance with the graded maximal exercise tests of the included studies (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Wilson et al., 1991;Marzo et al., 1992;Rudas et al., 1993;Kao et al., 1994;Doering et al., 1996;Geny et al., 1996;Notarius et al., 1998;Hayman et al., 2010;Peled et al., 2017;Nygaard et al., 2019;Nytrøen et al., 2019) and (2) the simulation duration should be sufficient to guarantee the stability of the model and therefore ensure the reliable determination of hemodynamic variables. Taking this into consideration, for resting condition, simulations of 300-s duration were performed. In the case of exercise simulations, the exercise intensities were increased in steps of 10% of maximum tolerable power until maximum exercise intensity was reached. As for resting condition, each exercise step was of 300-s duration. The simulated exercise intensity was adjusted so that peak CO complied with the pooled published data (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Wilson et al., 1991;Marzo et al., 1992;Rudas et al., 1993;Kao et al., 1994;Doering et al., 1996;Geny et al., 1996;Notarius et al., 1998;Hayman et al., 2010;Peled et al., 2017;Nygaard et al., 2019;Nytrøen et al., 2019). In HTxRs, exercise was terminated at 50% of the MTP of the control group, which is in accordance with the published clinical data (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Wilson et al., 1991;Marzo et al., 1992;Rudas et al., 1993;Kao et al., 1994;Doering et al., 1996;Geny et al., 1996;Notarius et al., 1998;Hayman et al., 2010;Peled et al., 2017;Nygaard et al., 2019;Nytrøen et al., 2019). For both rest and exercise conditions, mean values of hemodynamic variables were calculated from time intervals of 300-s duration for resting condition, and each exercise step. HR was derived by averaging the reciprocal time differences between the peaks of consecutive SA node action potentials, MAP as the average of aortic pressure, mean DBP, and mean SBP  from averaging minimum and maximum peak values of aortic pressures, respectively. CO was calculated as the average of left ventricular outflow, and stroke volume (SV) by averaging the differences between end-diastolic and end-systolic left ventricular volume. Systemic vascular resistance (SVR) was calculated as the average resistance of the systemic circulation compartment. Diastolic graft function was assessed through the end-diastolic volume index of the left ventricle (LVEDVI) and left ventricular end-diastolic pressure (LVEDP). The LVEDVI was computed by dividing the end-diastolic left ventricular volume by the mean body surface area (BSA) of the respective cohort. The BSA was calculated according to Gehan and George (1970).

Heart Rate
Simulations of resting condition predict baseline HR of 93.8 bpm (published data, 92.7 ± 10.3 bpm) for HTxRs. This equals an increase of about 25% compared to healthy baseline HR, which is consistent with the published data (Figure 3). Especially noticeable is the reduced chronotropic response in HTxRs. The model predicts a peak HR of only 133.9 bpm (published data, 129.3 ± 17.5 bpm), which represents a reduction of more than 20% compared to the healthy control group (Figure 3). We can clearly see the impaired chronotropic response in HTxRs, which is characterized by the notably lower slope of the curve compared to the control group ( Figure 4A).

Stroke Volume
In contrast to the healthy control group, simulations show a marked reduction in resting SV by 25% to only 55.4 mL (published data, 58.7 ± 15 mL) in HTxRs.
In response to exercise, simulations of HTxRs show strong SV augmentation, eventually reaching a maximum value of 95.7 mL (published data, 97.1 ± 19.7 mL), which is about 20% lower than in the healthy individuals (Figure 3).
The model predicts notably stronger SV augmentation in HTxRs than in age-and gender-matched healthy controls, although reaching subnormal values at their respective peak exercise ( Figure 4C).

Cardiac Output
Simulations of resting condition show normal basal CO of 5.2 L/min (published data, 5.9 ± 1.4 L/min) for HTxRs comparable to results for the healthy control group (Figure 3).
However, exercise simulations of HTxRs predict a peak CO of only 13 L/min (published data, 12.9 ± 3.1 L/min), which represents a drastic reduction of more than 40%, compared to healthy individuals (Figure 3). Despite the reduced peak CO, in response to low-to moderate-intensity exercise, the model predicts an almost identical increase of CO in both groups ( Figure 4B). Differences become clear only at moderate-to high-intensity exercise for which healthy persons can further increase their CO, eventually reaching nearly twice the values of HTxRs ( Figure 4B).

Left Ventricular End-Diastolic Volume
Consistent with the pooled published data, simulations of HTxRs at rest predict an LVEDVI of 51.9 mL/m 2 (published data, 49.2 ± 12.1 mL/m 2 ), which is about 20% lower than that of ageand gender-matched healthy individuals (Table 1).
FIGURE 4 | Comparison of exercise response between HTxRs and healthy controls, characterized by selected hemodynamic variables. The horizontal axis shows the relative exercise intensity with respect to the peak exercise intensity of the healthy control group (1,032 kpm/min). Note the substantial reduction of exercise capacity in HTxRs, reaching only about 50% (522 kpm/min) of the healthy control's peak exercise intensity. We can see a marked impairment of chronotropic response and reduced peak HR in HTxRs (A). Furthermore, the simulation results show no significant differences in augmentation in CO but substantially decreased peak CO in HTxRs (B). Resting SV is notably reduced in HTxRs and undergoes a rapid increase in response to exercise (C). MAP response is similar in both groups; however, healthy controls reach higher peak values due to greater exercise intensity (D). Despite the reduced absolute values, HTxRs show a markedly stronger augmentation of LVEDVI in response to exercise (E). Simulations predict a similar resting LVEDP for both the HTxRs, however, the augmentation of LVEDP is markedly stronger than in age-and gender-matched healthy individuals (F). HR, heart rate; CO, cardiac output; SV, stroke volume; MAP, mean arterial pressure; LVEDVI, left ventricular end-diastolic volume index; LVEDP, left ventricular end-diastolic pressure.
In contrast to the control group, simulations of exercise in HTxRs show a strong and immediate augmentation of LVEDVI, reaching 75.8 mL/m 2 at their respective peak exercise load, which is about 15% higher than that of healthy individuals at the same exercise intensity (Figure 4E).
Despite HTxRs show a stronger relative increase in LVEDVI, the model predicts notably higher total LVEDVI of 84.7 mL/m 2 at peak exercise in the control group.
However, only one publication (Kao et al., 1994) of the literature used for data pooling included values for LVEDVI at peak exercise. In this study, the peak exercise load was 390 and 825 kpm/min for HTxRs and the control group, respectively, which equals an exercise intensity of 37 and 71%, respectively, in our simulation. For these intensities, the model predicts an LVEDVI of 66.1 mL/m 2 (published data, 58 ± 11.6 mL/m 2 ) for HTxRs and 75.3 mL/m 2 (published data, 67.3 ± 15.1 mL/m 2 ) for age-and gender-matched healthy individuals, which is in good accordance with the values published in the literature (Kao et al., 1994).

Left Ventricular End-Diastolic Pressure
Simulations predict a similar resting LVEDP for both the HTxRs and the age-and gender-matched healthy control group with 2.8 and 3.6 mmHg, respectively ( Figure 4F).
In HTxRs, the augmentation of LVEDP is markedly stronger than in age-and gender-matched healthy individuals ( Figure 4F). HTxRs reach about 4.4 mmHg at peak exercise, which is about 20% higher than that of healthy individuals at the same exercise intensity. However, the LVEDP of 6.9 mmHg reached by healthy individuals at peak exercise is 30% higher than that of HTxRs at their respective peak exercise intensity.
None of the publications included for data pooling reported values of LVEDP at rest or exercise.

Model Validation
Except for the peak exercise MAP of healthy individuals, the simulation results are in good accordance with the pooled literature data, all being within the standard deviation of patient variability. Although the predicted peak exercise MAP is not in this range, that does not necessarily indicate a bad model representation of exercise hemodynamics. The small discrepancy between predicted and true MAP is consistent with the comparably lower SVR of the healthy study cohort (Braith et al., 1992) (Table 2). Overall, the simulation results are consistent with published clinical data (Pflugfelder et al., 1988;Hosenpud et al., 1989;Tamburino et al., 1989;Niset et al., 1991;Braith et al., 1992;Scott et al., 1995;Nytrøen et al., 2011;Cattermole et al., 2017), therefore verifying the validity of the presented model.
Simulations predict that the markedly increased resting HR in early HTxRs is exclusively due to the absence of vagal tone FIGURE 5 | Pressure-volume loops for a time period of 10 s obtained from simulations of healthy controls (A) and HTxRs (B), for rest (0%), medium (50%), and peak (100%) exercise intensity. Note that peak exercise in HTxRs corresponds to 50% of exercise intensity in healthy individuals.
resulting from parasympathetic cardiac denervation and that, supersensitivity to circulating catecholamines, is neglectable. However, normal CO is still maintained through a reduction of SV by about 25% compared to age-and gender-matched healthy individuals (Figure 3).
Simulations of maximum graded exercise do not show abnormally elevated pulmonary pressure in HTxRs at any intensity. HTxRs do reach a left atrial pressure (LAP) of 11 mmHg at peak exercise, which is similar to that of healthy controls, which attain a LAP of 11.6 mmHg at 50% exercise intensity, suggesting that reduced exercise capacity does not result from early dyspnea. However, for the sake of brevity, the results of the LAP response are not presented here.
The model predicts a strongly impeded chronotropic response in HTxRs, reaching only about 80% of the normal peak HR (Figure 3). The reduced HR response at low exercise intensities (up to 20%) results from a lack of vagal withdrawal due to complete vagal cardiac denervation. The inability to sufficiently increase the HR in response to exercise intensities greater than 30% is a consequence of the virtually absent sympathetic drive that is reduced by 95% due to sympathetic cardiac denervation. The remaining HR augmentation primarily relies on circulating catecholamines, which is also general consent in literature (Awad et al., 2016;Kobashigawa and Olymbios, 2017).
Ion-channel gating mechanisms of the SA node and βagonist clearing remained unchanged in HTxRs. Therefore, the model neglects the often assumed supersensitivity of the donor heart to circulating catecholamines (Awad et al., 2016;Kobashigawa and Olymbios, 2017). Nevertheless, model predictions for exercise response still are in excellent accordance with literature data, suggesting that supersensitivity may play a neglectable role in exercise response compensation in early HTxRs. A possible explanation is that supersensitivity typically tends to develop over time (Kobashigawa and Olymbios, 2017), thus being not present yet in the modeled population with a post-HTx time of 3.8 ± 2.9 months.
Moreover, simulations attain 20% less vasodilation in HTxRs (Figure 3), contributing to the limitation in CO augmentation. Impaired muscular vasodilation may be the result of prolonged periods of deconditioning pre-HTx (Kobashigawa and Olymbios, 2017) and changes in the vasopressor effect associated with cyclosporine treatment (Andreassen et al., 1998).
For low-to medium-intensity exercise, healthy individuals mainly rely on HR increase through vagal withdrawal to augment their CO (Figure 4), whereas enhanced ventricular contractility plays a minor role ( Figure 5A). However, for medium-to high-intensity exercise, healthy individuals further increase their CO primarily through SV augmentation, while onward HR augmentation plays a secondary role ( Figure 5A). The increase in SV can be attributed to raised contractility, characterized by the increased slope of the end-systolic pressure-volume relationship and increasing end-diastolic volume through the Frank-Starling mechanism ( Figure 5A).
Simulations show that despite the reduced exercise capacity, for low to moderate intensity, HTxRs can maintain normal CO, which is consistent with the findings in the literature (Kobashigawa and Olymbios, 2017). However, in contrast to healthy individuals, the exercise response in HTxRs is quite different. The model predicts that for any given exercise intensity, HTxRs strongly rely on the increase of preload through the Frank-Starling mechanism, enhancing their SV to compensate for chronotropic and inotropic incompetence ( Figure 5B). Nevertheless, an increase in HR through circulating catecholamines plays a non-neglectable, yet secondary role for CO augmentation (Figure 4).
Diastolic dysfunction is an often-reported condition in HTxRs manifesting in a leftward shift of the end-diastolic pressure-volume relationship (EDPVR); however, the underlying mechanisms are still not entirely understood. Possible reasons include cardiac denervation, heart graft remodeling, ischemia, and graft rejection (Rudas et al., 1990). In our model, we did not explicitly modify the EDPVR, thus not accounting for possible diastolic graft dysfunction. Consequently, the cardiac model represents a properly vascularized heart graft early after transplantation in which no considerable cardiac remodeling has taken place and that is not affected by graft rejection, which reflects the modeled population of early HTxRs (McLaughlin et al., 1978;Crisafulli et al., 1985;Kavanagh et al., 1988;Labovitz et al., 1989;Wilson et al., 1991;Marzo et al., 1992;Rudas et al., 1993;Kao et al., 1994;Doering et al., 1996;Geny et al., 1996;Notarius et al., 1998;Hayman et al., 2010;Peled et al., 2017;Nygaard et al., 2019;Nytrøen et al., 2019). Simulations to investigate the diastolic graft dysfunction and its influence on resting and especially on exercise hemodynamics should be the focus of future studies. Still, the model predicts evident acute changes in diastolic pressures and volumes in early HTxRs, showing markedly reduced resting and peak LVEDVI compared to the control group. The values for LVEDVI at rest and exercise are in good accordance with the literature (Kao et al., 1994;Nygaard et al., 2019), suggesting that despite the reduced absolute values, HTxRs show a markedly stronger augmentation of LVEDVI in response to exercise ( Figure 4E). Ultimately, the simulations predict similar LVEDP at rest in both groups, while the pressure increases notably stronger with exercise load in HTxRs than in the control group, allowing a greater increase in LVEDV.
The model was validated based on published clinical data of patients other than used for tuning of the model. For HTxRs, all predicted hemodynamic values are within the range of one standard deviation of the true values (Table 2), therefore suggesting a good model representation of resting and exercise hemodynamics in early HTxRs. Except for peak exercise MAP, the model predictions for the healthy group are also in good accordance with the published clinical data. The predicted peak exercise MAP is not in the range of one standard deviation of natural interpatient variability; however, the trend in the MAP response to exercise is the same as in the HTxR group. The small discrepancy in MAP is also consistent with the differences between predicted and true SVR ( Table 2). The sample size of the cohort used to validate the model of normal exercise hemodynamics is rather small (n = 10); therefore, it might be necessary to further validate the model based on hemodynamic data of a larger group.
Summarizing, according to the simulation results, in early HTxRs, the limited exercise response is primarily due to virtually absent sympathetic drive, which is reduced by 95% due to sympathetic cardiac denervation. The simulations highlight chronotropic incompetence following sympathetic cardiac denervation as the main limiting factor for exercise in HTxRs. Peripheral factors seem to play a secondary role, while pulmonary factors are negligible in the limitation of exercise capacity. The model shows that the main compensatory mechanism can be attributed to SV augmentation strongly relying on increased preload, while supersensitivity to circulating catecholamines is insignificant in early HTxRs.

Limitations of the Study
The model was validated with a small complementary set of clinical hemodynamic data of HTxRs and healthy individuals. However, the data contained only a limited number of patient measurements and hemodynamic parameters. Consequently, the model should undergo further validation, based on a more comprehensive dataset.
The often-observed diastolic dysfunction in HTxRs was not explicitly modeled. Thus, the cardiac model represents a properly vascularized, short-term post-HTx condition in which no remodeling has taken place yet. Simulations to investigate the diastolic graft dysfunction and its influence on resting and exercise hemodynamics were not part of this study but certainly interesting for future study, to investigate late HTxR hemodynamic response.

CONCLUSION
The present study provides a model-based perspective on the hypothesized origins of reduced exercise capacity and compensatory mechanisms for chronotropic and inotropic incompetence in early HTxRs.
Simulation results show an overall reduced exercise capacity in early HTxRs, which is primarily due to chronotropic incompetence, whereas peripheral factors play a secondary role. HTxRs can maintain normal CO for low-to medium-intensity exercise by compensation of chronotropic and inotropic incompetence through increased filling pressures.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
MH and FM contributed to the conception and design of the study and were involved in the interpretation of the study results. MH and DD conducted the literature review, contributed to the model development, and performed the simulations. DD performed the data pooling. MH wrote the manuscript and generated the illustrations. All authors contributed to manuscript revision, read, and approved the submitted version.