A Computational Study of the Effects of Tachycardia-Induced Remodeling on Calcium Wave Propagation in Rabbit Atrial Myocytes

In atrial cardiomyocytes without a well-developed T-tubule system, calcium diffuses from the periphery toward the center creating a centripetal wave pattern. During atrial fibrillation, rapid activation of atrial myocytes induces complex remodeling in diffusion properties that result in failure of calcium to propagate in a fully regenerative manner toward the center; a phenomenon termed “calcium silencing.” This has been observed in rabbit atrial myocytes after exposure to prolonged rapid pacing. Although experimental studies have pointed to possible mechanisms underlying calcium silencing, their individual effects and relative importance remain largely unknown. In this study we used computational modeling of the rabbit atrial cardiomyocyte to query the individual and combined effects of the proposed mechanisms leading to calcium silencing and abnormal calcium wave propagation. We employed a population of models obtained from a newly developed model of the rabbit atrial myocyte with spatial representation of intracellular calcium handling. We selected parameters in the model that represent experimentally observed cellular remodeling which have been implicated in calcium silencing, and scaled their values in the population to match experimental observations. In particular, we changed the maximum conductances of ICaL, INCX, and INaK, RyR open probability, RyR density, Serca2a density, and calcium buffering strength. We incorporated remodeling in a population of 16 models by independently varying parameters that reproduce experimentally observed cellular remodeling, and quantified the resulting alterations in calcium dynamics and wave propagation patterns. The results show a strong effect of ICaL in driving calcium silencing, with INCX, INaK, and RyR density also resulting in calcium silencing in some models. Calcium alternans was observed in some models where INCX and Serca2a density had been changed. Simultaneously incorporating changes in all remodeled parameters resulted in calcium silencing in all models, indicating the predominant role of decreasing ICaL in the population phenotype.

In atrial cardiomyocytes without a well-developed T-tubule system, calcium diffuses from the periphery toward the center creating a centripetal wave pattern. During atrial fibrillation, rapid activation of atrial myocytes induces complex remodeling in diffusion properties that result in failure of calcium to propagate in a fully regenerative manner toward the center; a phenomenon termed "calcium silencing." This has been observed in rabbit atrial myocytes after exposure to prolonged rapid pacing. Although experimental studies have pointed to possible mechanisms underlying calcium silencing, their individual effects and relative importance remain largely unknown. In this study we used computational modeling of the rabbit atrial cardiomyocyte to query the individual and combined effects of the proposed mechanisms leading to calcium silencing and abnormal calcium wave propagation. We employed a population of models obtained from a newly developed model of the rabbit atrial myocyte with spatial representation of intracellular calcium handling. We selected parameters in the model that represent experimentally observed cellular remodeling which have been implicated in calcium silencing, and scaled their values in the population to match experimental observations. In particular, we changed the maximum conductances of I CaL , I NCX , and I NaK , RyR open probability, RyR density, Serca2a density, and calcium buffering strength. We incorporated remodeling in a population of 16 models by independently varying parameters that reproduce experimentally observed cellular remodeling, and quantified the resulting alterations in calcium dynamics and wave propagation patterns. The results show a strong effect of I CaL in driving calcium silencing, with I NCX , I NaK , and RyR density also resulting in calcium silencing in some models. Calcium alternans was observed in some models where I NCX and Serca2a density had been changed. Simultaneously incorporating changes in all remodeled parameters resulted in calcium silencing in all models, indicating the predominant role of decreasing I CaL in the population phenotype.

INTRODUCTION
Models of cardiac electrophysiology are important tools for studying the mechanisms of heart disease and impaired cardiac cell function, and have been particularly relevant in the study of arrhythmic events such as atrial fibrillation (AF) (Trayanova, 2014;Grandi and Maleckar, 2016;Vagos et al., 2018). An important application of the models is interpretation and translation of knowledge gained from animal experiments performed on a variety of different species. Atrial cardiomyocytes have a less developed T-tubule system than their ventricular counterpart, and in many species this results in markedly different intracellular calcium dynamics. In both cell types calcium enters via L-type calcium channels, locally accumulates close to the membrane, and triggers calcium-induced calcium release from calcium release units (CRU) in the junctional sarcoplasmic reticulum (SR) into the cytosol. The complex T-tubular system of ventricular cells effectively ensures that every CRU is in close proximity with the membrane, and ensures that this process occurs simultaneously throughout the cell. In atrial cells without well-developed T-tubules, only a small fraction of the CRUs are close to the membrane and therefore activated by the initial L-type channel influx. Activation of the remaining CRUs occurs as a calcium wave which propagates from the membrane toward the center of the cell via a regenerative "fire-diffuse-fire" response. As one CRU is activated, calcium diffuses through the cytosol to cause a local concentration rise at the neighboring CRUs, which then reaches the threshold level for activation and releases more calcium and the process repeats resulting in a centripetally propagating calcium wave.
Given the heterogeneous nature of intracellular calcium signaling in atrial myocytes (Trafford et al., 2013), they are prone to developing calcium instabilities, such as local fluctuations in cytosolic calcium levels, which can translate into arrhythmogenic activity (Greiser et al., 2014;Voigt et al., 2014). Furthermore, sustained arrhythmic conditions such as AF can lead to severe electrophysiological remodeling of the atrial myocytes, including up-or down-regulation of numerous ion channels and other membrane transporters. These changes may affect the calcium cycling and propagation within the cell, and may potentially exacerbate repolarization instabilities. For instance, calcium alternans have been extensively observed and implicated as precursors of AF episodes in patients (Narayan et al., 2002).
It has also been shown in rat atrial myocytes that under basal conditions, calcium does not fully propagate toward the center of the cell, mainly due to extensive calcium buffering, and the diffusion barrier of the mitochondria and Serca2a (Bootman et al., 2011). It is believed that this lack of a calcium signal at the inner regions of the cell under normal conditions acts as a reserve to allow calcium release to be increased under conditions of higher contraction demand, such as during β-adrenergic stimulation.
Although atrial myocytes of rabbit and rat species share similar morphological traits, including size and the lack of a well-developed T-tubule system, it has been observed in line scans of rabbit atrial myocytes that calcium propagates in a fully regenerative manner (Greiser et al., 2014), unlike in rat atrial cells. However, the same study showed that replicating chronic AF conditions via rapid pacing of rabbit atrial myocytes resulted in a complete absence of calcium propagation inside the cell, a conditioned termed "calcium silencing." This behavior was suggested to be a cardio-protective mechanism of the cell to suppress arrhythmogenic after-depolarizations, and has been observed both in rabbit and in AF patient myocytes (Greiser et al., 2014). Additionally, several studies have pointed out the role of calcium buffering in regulating available free calcium in the cell, and ultimately being a promoter of calcium handling abnormalities (Greiser et al., 2014;Smith and Eisner, 2019). However, the exact mechanisms leading to calcium silencing, and their relative importance, remain unclear.
In the present paper we elucidate the underlying mechanisms of calcium silencing by comparing computational model results with experimental observations. We have previously developed a model of the rabbit atrial myocyte with spatial calcium description (Vagos et al., 2020), which was parameterized to match experimentally observed rabbit-specific electrophysiology using a population of models approach. The parameterization resulted in 16 different models of the rabbit atrial cell, each with a unique combination of ion channel maximum conductances, but all closely matching experimentally observed electrophysiological characteristics. All models showed a rabbit-like action potential (AP) and calcium transient (CaT), as well as fully propagating calcium waves. The purpose of the present study is to assess the role of different cellular parameters that are believed to change as a result of rapid atrial pacing (RAP) induced remodeling. By changing their values in the model in one-at-a-time fashion, and applying these changes to all of the 16 models obtained from our previous study (Vagos et al., 2020). We changed the maximum conductances of I CaL , I NCX , and I NaK , RyR open probability parameters, RyR density, Serca2a density, as well as calcium buffering strength. We chose these specific parameters since all of them were shown to be altered in rabbit atrial myocytes as a consequence of RAP remodeling (Greiser et al., 2014). Furthermore, we scaled these parameters in the model according to experimental data reported in Greiser et al. (2014). Therefore, the changes introduced in the model to incorporate remodeling effects are based on direct measurements rather than indirect estimates from AP or CaT observations. Furthermore, it has been extensively shown that cardiomyocytes have a large degree of variability in ion channel expression levels and electrophysiological properties. Therefore, the advantage of changing the parameters in 16 different models is that it allows to incorporate electrophysiological variability in the simulations, providing a more comprehensive assessment of the effects of each individual parameter. Using a population of models rather than a single "representative" model, we obtain a more robust and insightful analysis of the model behavior and associated uncertainty. By testing the effects of parameter changes in different model instances it is possible to obtain a more generalized response of the model and to identify possible drivers of the observed responses.

Control Model Population
We used the previously published model of the rabbit atrial cardiomyocyte with spatial description of the calcium handling system, which allows for the simulation of intracellular calcium wave propagation. A complete description of the model development and structure is provided in Vagos et al. (2020). In brief, the model was based on the rabbit atrial myocyte model by Aslanidi et al. (2009), and the human atrial myocyte model by Voigt et al. (2014). It is composed of discrete units, each containing ionic concentrations and cellular structures specific of its location in the cell. Conceptually, the model is based on segmenting the cardiomyocyte into transversal slices, which are further sub-divided into domains in the transversal direction, as illustrated in Figure 1. As described in detail in Vagos et al. (2020), the model distinguishes between the domains close to the membrane and the interior domains. The inner domains contain the cytosol, sarcoplasmic reticulum (SR), a sub-SR space (SRS), and calcium release units (CRU), while the membrane domains also include the sub-sarcolemmal (SL) space, junctional space (j), and associated sarcolemmal currents. The membrane currents include calcium currents I CaL and I CaT ; the fast sodium current I Na ; repolarizing K + currents I to1 , I Kr , I Ks , and I K1 , three background currents I Cab , I Nab , and I Clb , as well as the Na + -Ca 2+ exchanger I NCX , the Na + -K + pump I NaK , and the plasmalemmal Ca 2+ pump I CaP . The formulation of the I NCX current was adopted from Voigt et al. (2014), while all other ionic current formulations were taken from Aslanidi et al. (2009), and reparameterized in Vagos et al. (2020). Flux of Ca 2+ between neighboring domains is described by a linear diffusion model adopted from Voigt et al. (2014). The model equations are specified in detail in Vagos et al. (2020), and the source code is available 1 .
The calibrated population of models developed in Vagos et al. (2020) consisted of 16 models with rabbit-like electrophysiological parameters: (1) fully regenerative calcium wave propagation from periphery to center of the cell; (2) absence of alternans or afterdepolarizations; and (3) AP and whole cell CaT morphologies within the physiological ranges experimentally observed in rabbit atrial myocytes. The models differed only in the values of maximum conductances of the ionic currents, thus capturing the natural variability observed in atrial cells. In the following we will refer to these 16 models as the "control population, " and its main characteristics are illustrated in Figure 2. Figure 2A shows APs for all 16 models, Figure 2B shows the CaTs from the membrane domains (CaT m ) and the innermost domains (CaT c ), while Figures 2C,D show an exemplary AP, CaTs, and calcium wave from one of the models. We refer to Vagos et al. (2020)

Incorporating RAP Remodeling
An experimental study has shown that rabbit atrial myocytes subjected to rapid atrial pacing (RAP), at 10 Hz for 5 days, developed an array of cellular alterations at the structural and biochemical levels (Greiser et al., 2014). Most notably, observations included decreased (1) I CaL and (2) I NaK currents, (3) increased I NCX , (4) reduced Serca2a density, (5) increased ryanodine receptor (RyR) open probability, (6) decreased RyR cluster density, and (7) increased calcium buffering strength. The membrane currents and Serca2a density were modified by g the maximum fluxes, to increase the RyR open probability we scaled parameters RyR_P[0], RyR_P[1], and RyR_P[5], as suggested in Voigt et al. (2014), decreased RyR density was incorporated by scaling RyR_P[11] and the NRyRs parameter, while buffering was increased by a scaling of the parameter Buff_factor. The parameter changes are specified in Table 1, and we refer to Vagos et al. (2020) and Voigt et al. (2014) for a specification of the model equations and parameters. These cellular alterations have been proposed as being a cardioprotective mechanism against arrhythmogenic calcium activity (Greiser, 2017). We introduced these alterations in the control model population defined above to assess the effects of each alteration on the action potential (AP) and calcium wave propagation dynamics, and also the effect of all the changes combined. The remodeling was incorporated by scaling the baseline values of the model parameters that modulate the affected cellular mechanisms to match experimental measurements. Each of the seven alterations listed above were implemented individually to the 16 models of the control population, and we also created one population that incorporated all the listed modifications collectively, resulting in eight remodeled populations with 16 models in each. The parameter changes applied for each of the remodeled populations are specified in Table 1. By applying the remodeling to a calibrated control population rather than a single representative model, we were able to capture more of the natural variability and parameter uncertainty and provide a more complete picture of the potential effect on calcium dynamics.

Visualization and Analysis of Calcium Waves
Calcium wave propagation in the models was visualized in spatio-temporal plots of CaTs from each cell domain over time, similar to line scan images of individual myocytes. We also characterized individual remodeling effects by studying the occurrence of CaT amplitude alternans and calcium silencing. Calcium alternans arise from instabilities in the calcium handling system (Weiss et al., 2006;Gaeta and Christini, 2012;Xie et al., 2014), and have also been clinically associated with atrial arrhythmia (Franz et al., 2012). We defined alternans as beat-to-beat differences in CaT amplitude larger than 5%, while calcium silencing was defined as a ratio of central to membrane CaT amplitude (CaT c amplitude/CaT m amplitude) ≤ 0.10 for all beats. Based on the presence of alternans and calcium silencing, the individual models were classified into five different categories:
Finally, we also report APD90 and resting memembrane potential (RMP) for the remodeled populations, since APD90 shortening and RMP depolarization have been linked to arrhythmogenic activity in cardiomyocytes (Bosch et al., 1999;

RESULTS
In this section we list the effects of the remodeling alterations defined in Table 2, including both the individual mechanisms and the combined effect of all mechanisms. The results are summarized in Figure 3 and Table 3. Figure 3 includes APs and CaTs for the eight model populations defined in section 2.2. The individual changes observed in each population are described in more detail below. Table 3 classifies each population according to the characteristics described in section 2.3, while Figure 4 shows spatio-temporal plots of selected models to illustrate the characteristics of each category. For each row in Figure 4, the left panel shows APs and CaTs from the subsarcolemmal (CaT m ) and central domain (CaT c ), while the right panel shows spatio-temporal linescan plots illustrating the calcium propagation. Figure 4A show normal calcium propagation (category 1), Figures 4B,C show examples of calcium silencing (category 3), Figure 4D shows alternans (category 4), Figure 4E displays calcium silencing and alternans (cat 5), and Figure 4F shows a pattern we have characterized as unphysiological (category 5, "other"). The example in Figure 4A is from the control population, Figure 4B is from population 1 (reduced G CaL ), Figure 4C from population 3 (reduced I max NaK ), Figure 4D from population 4 (reduced Serca2a density), Figure 4E from population 2 (increased I max NCX ), and Figure 4F is from population 5 (increased RyR open probability). The number of models falling into each category is quantified in Table 3. The table shows the percentage of models from each of the eight populations that fall into each category of calcium propagation. Figure 3A shows the population after reduction of G CaL to 40% in the population, which resulted in a consistent response of calcium silencing across all models. An example of a calcium wave from this population is shown in Figure 4B. Additionally, CaT m showed a considerably reduced amplitude compared to the control population, as would be expected since the amplitude of the sarcolemmal CaT is directly linked to the magnitude of I CaL . We also observed a significantly shortening of APD90 (see Table 2), which is a well-documented effect of reduced I CaL in myocytes.

Up-Regulation of I NCX
The effects of increasing I max NCX in the population resulted in a variety of calcium propagation patterns, including calcium alternans, silencing, and alternans and silencing (see Figure 3B). Additionally, this population showed large variations in APD90, with a mean APD90 larger than in the control population, as well as large variations in CaT amplitudes and diastolic calcium levels. Furthermore, one of the models failed to repolarize, as can be seen from the ripples in Figure 3B, left panel, with corresponding rippled CaTs (not shown). This model corresponded to a [Na + ] i level of 41 mM, and [Ca 2+ ] i of 5.8 µM. The perturbations in AP and CaTs observed in this population are not surprising given the known role of I NCX in regulating both [Ca 2+ ] i and repolarization (Weber et al., 2001).

Decreasing I max
NaK resulted in about 60% of the models developing calcium silencing, while the remainder did not change significantly. This population showed significantly elevated [Na + ] i levels (22 ± 4 mM vs. 9.5±2 mM, (two-sample t-test, α = 0.05), as expected given the decreased outflow of Na + from the cell, which shifted the activity of I NCX . As a consequence of reducing I max NaK , the forward mode (I NCX,fwd ) and reverse mode (I NCX,rev ) components of I NCX were significantly smaller and larger, respectively, in the I max NaK population as compared to the control (two-sample t-test, α = 0.05). This resulted in an accumulation of calcium in the cytosol and in the SR. The RMP was also significantly more depolarized in this population.
The CaT m amplitude and diastolic calcium levels also varied considerably in this population (see Figure 3C and Table 2). There were no observable differences in the maximum conductances of the models that showed normal propagation versus calcium silencing, but found that models showing calcium silencing mostly corresponded to higher values of both [Na + ] i and [Ca 2+ ] SR .

Decreased Serca2a Density
Decreasing Serca2a density to 70% of baseline value resulted in largely unchanged APs, and an overall small effect on CaT m and  CaT c , with only about 20% of the models showing alternans and the others being mostly unaltered.

Increased RyR Open Probability
Increasing RyR open probability resulted in about 40% of the models developing calcium silencing. We did not observed any apparent associations between the occurrence of alternans in this population, and the values of the altered parameters, or of [Na + ] i and [Ca 2+ ] SR .

Decreased RyR Density
The effect of decreasing RyR density largely resulted in a variety of unphysiological calcium propagation patterns. One such wave pattern is shown in Figure 4F. This result indicates that the calcium system in this model is quite sensitive to such large changes of the two target parameters, NRyRs and RyR_P[11].

Increased calcium Buffering Strength
Increasing the calcium buffering strength in the population did not result in any form of abnormal calcium wave propagation, with all models showing a normal calcium propagation pattern. This is in contrast with experimental observations from Greiser et al. where the authors reported that an increased calcium buffering strength in the cell was likely to be correlated with the observed calcium silencing in the remodeled cells.

Full RAP Remodeling
Finally, we simultaneously incorporated all seven RAP remodeling mechanisms considered above. This change, incorporating all 10 parameter changes listed in Table 1, resulted in a population with significantly shorter APD90, reduced CaT m amplitude, and calcium silencing. This result is in good agreement with experimental observations in Greiser et al. (2014), where the authors observed calcium silencing in line scans of RAP remodeled cells, without the occurrence of alternans or other pathological calcium propagation behaviors.

DISCUSSION
We used a previously developed population of 16 models representing normal rabbit atrial electrophysiology, and incorporated known remodeling effects from rapid atrial pacing. Seven new populations were built by changing parameters to reflect isolated physiological alterations, such as altered expression levels of an individual ionic current. Finally, one population incorporated the changes in all parameters simultaneously to reproduced full RAP-induced remodeling, as has been observed experimentally (Greiser et al., 2014). The approach used here to quantify the effects of changing model parameters on model behavior differs from the more traditional approach where a single model is used to represent a certain region of the heart, animal species, or genotype. In contrast, we opted to use a small population of 16 models, obtained by calibration of a population of 3,000 models to experimentally measured rabbit electrophysiological parameters. An advantage of the population approach over a single model approach is that it allows to incorporate uncertainty and natural variability observed in the atrial myocytes, and may provide a more robust analysis of the effects of electrophysiological remodeling.
In our simulations each remodeled population showed different phenotypes, with some parameter changes leading to consistent changes across the entire population, while others resulting in a myriad of different phenotypes. The overall effects of individual parameter changes on calcium wave patterns are compiled in Table 3, while Table 2 summarizes the effects on APD90, RMP, and CaT in all model populations. Mean CaT m amplitude was significantly reduced in the populations with altered G CaL , I max NCX , I max NaK , and RyR_P [11], while CaT c amplitude was reduced in all populations, except the one with altered calcium buffering. These amplitude changes resulted in an overall reduced CaT ratio in the populations with altered G CaL , I max NCX , I max NaK , Serca2a density, RyR Po, and RyR density.
In our simulations calcium silencing was mainly driven by reduction of G CaL and I max NaK , when these parameters are changed individually. The similarities between the G CaL and full remodeling populations indicate that the 60% reduction of G CaL was the dominant effector of changes in AP and CaT morphologies. Furthermore, despite significant differences in APD90 and RMP among the various remodeled populations, with one model showing an afterdepolarization and another showing failed repolarization, APs were mostly consistent across the populations. Greater differences were observed in CaT m and CaT c morphologies, especially in the G CaL , I max NCX , RyR Po, RyR density, and full remodeling populations. Furthermore, modification of G CaL , Serca2a density, Buff_factor, and full remodeling resulted in consistent alterations in calcium wave propagation patterns, while changes in the other parameters entailed a more complex response, despite considerable differences in the values of their maximum ionic conductances. In contrast, changes in the other parameters which resulted in a more varied set of responses indicate that those parameters play a more complex role in the calcium handling system, where the values of maximum ionic conductances may also have had an impact.
Experimental observations from Greiser et al. (2014) suggested that the main driver of calcium silencing in rabbit atrial cells could be increased calcium buffering strength, whereas FIGURE 4 | Exemplary calcium waves from the RAP remodeled populations showing calcium propagation patterns from the 5 wave categories. The left column shows AP and calcium transients, while the right column shows the corresponding spatio-temporal calcium wave plots. (A) corresponds to category 1 (normal), (B,C) to category 3, (D) to category 2, (E) to category 4, and (F) to category 5.
Frontiers in Physiology | www.frontiersin.org the reduction in I CaL was unlikely to play a significant role. Our computer simulations of the full remodeling case, which aims to emulate the array of cellular alterations measured experimentally, are consistent with the experimental observations of calcium silencing following RAP remodeling. However, our results indicate that reduced I CaL , rather than increased calcium buffering strength, may be the main effector of calcium silencing. calcium alternans were only observed in two of the remodeled populations (I max NCX and Serca2a density), and in a fairly low percentage of the models. This low occurrence of alternans suggests that they are the result of more drastic changes in cellular mechanisms, which were not captured in the remodeled populations constructed here. Overall, the results presented here provide insight into the possible mechanisms of calcium silencing in rabbit atrial myocytes.
As with any modeling study, the approach developed here has a number of limitations and assumptions. We note that, although the approach proposed here, whereby variability in maximum ionic conductances is taken into account, can improve robustness and reliability of simulation results, the control population used in this study is only an abstraction for a real population of myocytes obtained from different specimens and varying regions of the heart. The control population was produced by experimentally calibrating a population of 3,000 models against rabbit AP and CaT characteristics (Vagos et al., 2020), and thus the variation in maximum conductances introduced in the population cannot be verified experimentally as of this point due to the lack of data. A more holistic approach would be to take into account experimentally observed variability in ionic channel expression levels, for instance, and to calibrate the simulated population accordingly. The populations could also be expanded via re-sampling of the maximum conductances, such as introducing perturbations around the baseline values, and recalibrating.
While the RAP populations were based on previously published changes in ionic currents and calcium dynamics during RAP and AF, there is considerable variability in the reported changes, see for instance (Greiser et al., 2011). Our modifications are mainly based on experimental data from Greiser et al. (2014), and while all the modifications are supported by experimental data they are by no means the only possible alternative. A natural extension of the present work would be to extend the population of models approach to consider a wider range of previously reported electrophysiological remodeling, to capture more of the underlying uncertainty and potentially gain insight into plausible mechanisms underlying the observations. Another very natural extension would be to incorporate variability in the RyR parameters as well as J max up and Buff_factor, thereby producing a population that more reliably captures the full spectrum of calcium wave propagation patterns resulting from RAP remodeling. Such an extension could also be complemented by a more detailed quantification and characterization of the calcium propagation properties, providing more detailed insight than the broad categorization performed here.

CONCLUSIONS
In this study we have proposed an approach for analyzing the effects of rapid atrial pacing remodeling in a population of healthy rabbit atrial myocytes. The results indicate that reduced I CaL is the main determinant of the occurrence of calcium silencing, although decreased I max NCX , increased I max NaK , increased RyR open probability, and reduced RyR density also resulting in some cases of calcium silencing. Reduction in I NCX , and in Serca2a density gave rise to calcium alternans in a small percentage of instances, while changes in RyR density largely resulted in unphysiological calcium transients. The methodology proposed in this paper provides a robust framework for assessing effects of cellular remodeling under conditions of uncertainty and natural variability, and offers additional insight compared with studies based on a single representative model.

DATA AVAILABILITY STATEMENT
Parameter sets and simulation data presented in this study can be found at: https://github.com/marciavagos/rabbit_model_RAP_ data.git.

AUTHOR CONTRIBUTIONS
MV contributed to conceiving the study, implemented the model, performed all numerical simulations and data analyses, wrote the first draft of all sections and revised, and edited all sections. JS and HA contributed to conceiving the study, contributed to the analysis and model development, and revised and edited all sections of the paper. US and JH contributed to conceiving the study and revised and edited all sections of the paper. All authors contributed to the article and approved the submitted version.

FUNDING
This project has received funding from the European Union's Horizon 2020 MSCA ITN under grant agreement No. 675351.