Multiscale Modeling of Dyadic Structure-Function Relation in Ventricular Cardiac Myocytes.

Cardiovascular disease is often related to defects of subcellular components in cardiac myocytes, specifically in the dyadic cleft, which include changes in cleft geometry and channel placement. Modeling of these pathological changes requires both spatially resolved cleft as well as whole cell level descriptions. We use a multiscale model to create dyadic structure-function relationships to explore the impact of molecular changes on whole cell electrophysiology and calcium cycling. This multiscale model incorporates stochastic simulation of individual L-type calcium channels and ryanodine receptor channels, spatially detailed concentration dynamics in dyadic clefts, rabbit membrane potential dynamics, and a system of partial differential equations for myoplasmic and lumenal free Ca2+ and Ca2+-binding molecules in the bulk of the cell. We found action potential duration, systolic, and diastolic [Ca2+] to respond most sensitively to changes in L-type calcium channel current. The ryanodine receptor channel cluster structure inside dyadic clefts was found to affect all biomarkers investigated. The shape of clusters observed in experiments by Jayasinghe et al. and channel density within the cluster (characterized by mean occupancy) showed the strongest correlation to the effects on biomarkers.

Significance

Diseases such as myocardial infarction, aortic stenosis, tachycardia, hypertension, chronic ischemia, and atrial fibrillation have been related to changes inside the dyadic cleft, which is a subvolume of cardiac myocytes of ∼10−17 l (typical cell volume 10−11 l). However, exploration of the relation between subdyadic structures and disease is difficult because such microscopic structures in cells are in many cases not amenable to experimental manipulation, or experiments addressing them might not allow for simultaneous observation of cellular responses. Multiscale mathematical models can explore the relation between microscopic structures and cellular response. We show by mathematical modeling that the geometric properties of ryanodine receptor channel clusters within dyadic clefts affect cellular responses.

Introduction

The functioning of the heart is based on the precisely controlled contraction of its cardiac myocytes coordinated across the muscle by waves of membrane potential depolarizations (action potentials (APs)) emanating from the sinoatrial node. On the level of an individual cardiac myocyte, L-type Ca2+ channels (LCCs) open during an AP and trigger the release of Ca2+ from the sarcoplasmic reticulum (SR), which is the main intracellular Ca2+ storage compartment. The entailing global Ca2+ increase causes the binding of molecular motors to actin filaments in the sarcomeres and initiates contraction.

The SR forms a network of tubes that extend throughout the interior of the cell and can be divided into two main components known as junctional SR (jSR) and network SR (nSR). The cardiac myocyte is penetrated by a network of transverse tubules, which are plasma membrane invaginations that approach the jSR and thereby form small cellular subvolumes (see Fig. 1). These subvolumes, which have a height of 10–15 nm, are called dyadic clefts. Action potential gated LCC opening leads to calcium-induced calcium release (CICR) through ryanodine receptor channels (RyRs) in the jSR membrane. Ca2+ induces its own release because the opening probability of the RyRs increases with the local Ca2+ concentration at the channel. CICR is facilitated by the co-localization of LCCs in the transverse tubules membrane on the one side of a dyadic cleft and RyRs in the jSR membrane on the other side. The RyRs together with the LCCs and the associated jSR structure comprise the calcium release unit (CRU). Release from CRUs provides the Ca2+ triggering contraction of the sarcomeres. CRUs are concentrated in the structures connecting sarcomeres—the z-disks. The CRUs are arranged within z-disks with distances of less than 1 μm. The z-disks form a regular stack with a spacing of ∼2 μm.

Figure 1

(A) Sketch of cellular level organization and transverse tubular structure in the cardiac myocyte. (B) Molecular level arrangement of LCCs and RyRs in a single dyadic cleft is shown. To see this figure in color, go online.

CRUs behave stochastically because they contain a small number of ion channels. Because of CICR, they are excitable and can form sparks, which are the elementary events of Ca2+ release. Cooperation of several CRUs via CICR may generate unwanted Ca2+ waves. These processes occur on different time and length scales (see Fig. S1). The [Ca2+] changes inside the dyadic cleft happen within a few milliseconds; SR dynamics, on the other hand, act on a timescale of up to tens of seconds. We have spatial scales ranging from tens of nanometers in the dyadic cleft up to 100 μm in cell size. To account for these temporal and spatial scales, multiscale models with spatially distributed Ca2+ release sites have been developed (1, 2, 3, 4, 5, 6). The model used here simulates the behavior of individual RyR and L-type calcium channels as well as the concentration gradients inside clefts. On cell level, we simulate the membrane potential and concentration dynamics. We do not use the approximation by spatial compartments for the bulk concentration dynamics but simulate the corresponding partial differential equations with the numerically required spatial resolution (1, 6).

Common challenges with detailed multiscale modeling are the parameterization of the model by reproducing the values of a set of measured biomarkers (biomarkers are measurable properties of cells characterizing cell behavior) and the quantification of unknown values of parameters (7, 8, 9, 10) because of the large number of simulations required for these purposes. Methods to quantify the relationship of variability and uncertainty of model inputs (parameter values) to outputs (simulated biomarker values), which are based on the construction of a response surface, have been recently suggested (7, 8, 11, 12) and may spare many simulation runs. Exploiting the multiscale abilities of the model, we use biomarkers on different time and length scales for APs (seven biomarkers) and Ca2+ sparks (four biomarkers) to adapt our model to rabbit experimental results. Our systematic quantification of the relation between parameter and biomarker values uses an approximation by sums of polynomials (polynomial chaos expansion (12, 13, 14)). Sensitivity analysis identifies the parameters dominating the control of biomarker values—both in the mathematical sense of it and as the model’s suggestion for most efficient control of the cell state (e.g., by post-translational modifications).

Diseases such as myocardial infarction, aortic stenosis, tachycardia, hypertension, and chronic ischemia are frequently related to changes in the dyadic cleft (15, 16), which motivated several studies in recent years focusing on the details of the placement of RyR channels inside it. It turned out not to be on a square lattice as assumed before but to be less regular with respect to size and geometrical properties (16, 17, 18, 19, 20, 21, 22, 23, 24). The geometrical analysis revealed that channel positions in a cluster have random components and that cluster area is elongated in one direction rather than quadratic or circular (17, 18, 19, 20, 24, 25). Although the cluster size heterogeneity has been related to spark probability (23), the functional consequences of the geometrical properties of RyR clusters are not so obvious yet. Modeling can investigate them only if intradyadic gradients are taken into account as our approach does and other studies on the CRU level did (22, 25, 26). We address the functional consequences of channel placement on CRU and cell level and compare them between the regular arrangement and configurations with increasing irregularity. To that end, we choose the rules of channel placement provided by Jayasinghe et al. (17) to generate cluster geometries and channel locations similar to experimental observation. We put CRUs with these cluster geometries into a ventricular cell model to study the relation between cluster structure and cell function.

Materials And Methods

Mathematical model and methods

The mathematical model comprises whole cell dynamics as well as local molecular events (see Figs. 1 and S1). On the finest level, individual channels are represented as continuous time Markov chains, coupled by local gradients inside the dyadic space. Cell wide diffusion of [Ca2+] and its buffers is modeled by partial differential equations, which also include the fluxes generated by SR/endoplasmic reticulum Ca2+-ATPase (SERCA), NCX (Na+/Ca2+ exchanger), and the CRUs. Spatially averaged variables comprise membrane potential, [Na+], and [K+] which are generally assumed to not exhibit strong gradients on the subcellular level. A detailed description of the model can be found in (1, 6) and Supporting Materials and Methods.

The finite element simulation toolbox DUNE has been used to solve the model equations (27, 28). A complete description of the numerical approach is given in (1, 6) and a short overview in Supplementary Materials and Methods.

Channel placement model

Mathematical models describing cellular calcium dynamics in cardiac myocytes generally omit a spatially resolved description of the single CRUs depicting them as point sources and/or neglect the internal structure of channel arrays inside dyadic clefts by assuming spatially homogeneous Ca2+ concentration (4, 29, 30, 31). However, we take them into account as required by our investigation on RyR placement (1, 6, 32, 33). To determine the RyR arrangement in the cleft, we use the placement algorithm suggested by Jayasinghe et al. (17), which provides channel locations closely resembling their experimental data. This placement algorithm determines channel locations as the sequence of positions in a two dimensional random walk with as many steps as channels in the dyad. The first RyR is placed in the center of the dyadic cleft. The position of the second RyR is found by a step in a random direction and with random length. The step length is drawn from a normal distribution (mean μRyR = 40.1 nm, SD σRyR = 7.4 nm from (17)) with a cutoff accounting for the channel molecule diameter of 30 nm (34). The angle defining the direction is drawn from a uniform distribution in [0,2π]. Subsequent steps to channel positions obey the same rules plus the additional requirement to steer clear of existing channel molecules (excluded volume).

The first LCC channel is positioned at the center of the RyR cluster. In the following, LCCs are placed on a regular grid as in (1), again with a minimal distance of 30 nm from any other channel. Two examples of channel locations in a dyadic cleft generated this way are shown in Fig. 2 C.

Figure 2

(A) L-type Ca2+ channel (LCC) state scheme with the open state marked in green (36). (B) Ryanodine receptor (RyR) state scheme with the open state marked in green (26); the open rate kopen depends on dyadic as well as on jSR [Ca2+] as described in more detail in Supporting Materials and Methods. (C) Shown are two examples of placements of LCCs (blue circles) and RyRs (red diamonds) in single dyadic clefts with radii r. To see this figure in color, go online.

Sensitivity analysis and construction of a response surface

We generated a population of models by varying five crucial model parameters using Latin hypercube sampling (35). The hypercube was formed by the axes in the parameter space representing kplus, kclose, gRyR, gLCC, and VP,max. All five parameters were varied by a factor of 10 (see Table 1). The choice of the parameter ranges was based on values in the literature (1, 36, 37). The literature values may depart slightly from the chosen ranges because of different pacing cycle lengths. The pacing cycle length used in our simulations is 350 ms. The model simulations have identical initial conditions except for the stochasticity of the geometric channel arrangement and selected model parameters. All samples were run for the same simulation time. The resulting set of simulation results was analyzed by Bayesian linear regression (38) and polynomial chaos expansion (13) to obtain local and global parameter dependencies.

Table 1

Parameters, Sampling Ranges, and Suggested Parameter Set for Spark and Action Potential Simulations

Parameter	Description	Sampling Range	Suggested Value	Accepted Range
kplus [ms−1μM−η]	RyR opening rate	5.0 × 10−5–5.0 × 10−4	1.5 × 10−4	(1.1–2.3) × 10−4
kclose [ms−1]	RyR closing rate	0.1–1.0	0.5	0.28–0.55
gRyR [μm3 s−1]	RyR Ca2+ permeability	3.0 × 10−4–3.0 × 10−3	7.5 × 10−4	(5.9–8.4) × 10−4
gLCC [μm3 s−1]	LCC Ca2+ permeability	4.5 × 10−4–4.5 × 10−3	3.2 × 10−3	(1.7–3.4) × 10−3
VP,max [μM ms−1]	maximal SERCA uptake rate	0.15–1.5	0.55	0.08–0.71

The output stays in the literature range of all biomarker values, if the corresponding parameter is varied within the range given in the fifth column (accepted range), whereas other parameter values are kept at the value in the fourth column (suggested value).

Bayesian linear regression was used to obtain an estimate of local parameter sensitivity coefficients. The model sensitivity s from the linear fit was computed from the following:where S is the slope (along parameter X) from multiple linear regression in a specified neighborhood of a reference parameter set (marked in Figs. 5 and 7 as a red cross). Here Xref denotes the corresponding parameter value, and Yref is the corresponding reference value of the biomarker. We can read off the local strength, direction, and uncertainty of output change with respect to the variation of a selected parameter from these sensitivity coefficients. Because the mathematical model is stochastic, the output cannot be predicted with absolute certainty but with some probability only. We therefore used Bayesian linear regression to obtain a proper quantification of the uncertainty in the predictions. The algorithms were implemented using the Python library Edward (39) and TensorFlow (40) toolboxes. Although these sensitivities are based on a linear regression, we complemented our investigations by calculating Sobol coefficients (see Fig. S3), which serve as a measure of global sensitivities for nonlinear models (41).

Figure 5

(A–C) Contour plots for the mean values of APD90, peak systolic , and mean diastolic in dependence on LCC permeability gLCC and maximal SERCA uptake VP,max with the values for kplus, kclose, and gRyR fixed to the values in Table 1 (fourth column). The iso lines for the upper and lower parameter values of the literature ranges are color coded in green and orange, respectively. (D) Shown are the contours limiting the literature value ranges for all three biomarkers in a single plot. The white area outlines the parameter region for gLCC and VP,max within which all three AP biomarkers are within the literature ranges. The red mark indicates the parameter set listed in Table 1 (fourth column), meeting also the spark biomarker requirements. To see this figure in color, go online.

Figure 7

(A–C) Contour plots for the mean values of FDHM, spark rate, and peak inferred from Fluo-4 bound Ca2+ (Eq. 3) in dependence on the RyR permeability gRyR and opening rate kplus with the values for gLCC, kclose, and VP,max fixed to Table 1 (fourth column). The iso lines for the upper and lower parameter values of the literature ranges are color coded in green and orange, respectively. (D) Shown are the contours limiting the literature value ranges for all three biomarkers in a single plot. The white area outlines the parameter region for gRyR and kplus within which all three spark biomarkers are within the literature ranges. The red mark indicates the parameter set listed in Table 1 (fourth column), meeting also the AP biomarker requirements. To see this figure in color, go online.

To quantify the effect parameters have on the biomarkers, we used an approximation method known as polynomial chaos expansion. For the polynomial chaos expansion, an orthonormal basis of polynomials was generated by using a Python-specific library called Chaospy (42). We assumed uniform distribution of the input parameters and therefore used Legendre polynomials for the regression fit. Using a polynomial degree p and a number of parameters d, the number of polynomial coefficients, which have to be determined, can be calculated from:

The required number of data points usually exceeds the number of coefficients by at least a factor of 2–3 to prevent overfitting (43). To obtain an optimal regression and polynomial degree, we quantified the commonly used least-squares fit error and the cross-validation error as explained in Supporting Materials and Methods.

Results

The simulated time course of the Ca2+ concentration inside a dyadic cleft during an AP is illustrated in Fig. 3 and Video S1. Gradients comprise three orders of magnitude (0.1–150 μM) upon the opening of the first channel. The concentration outside the cleft space changes quickly, such that we observe gradients from ∼150 μM at the boundary of the cleft to ∼300 μM at open channels later during the event. Hence, the opening rate of RyRs close to open channels is initially six orders of magnitude and, later into the spark, is ∼4 times faster than the rate of channels further away (see Eq. S11; Table S2). The simulations illustrate the strong impact of gradients on the transition from one open LCC (or RyR, quark) to a spark and the quantitative effect on the calcium transient. This applies throughout an AP, as Video S1 shows.

Figure 3

Multiscale simulation. Shown is a snapshot of the concentration profile in a dyadic cleft and the jSR concentration time course of this CRU. A snapshot of isoconcentration surfaces (0.6 μM, 2.6 μM) of the cytosolic concentration [Ca2+]i of the whole z-disk (15 μm × 15 μm × 2 μm) with 320 dyadic clefts and the time course of the membrane potential, average [Ca2+]i, and average nSR [Ca2+] are shown on the right-hand side. The red line indicates the time point at which the snapshots were captured. A corresponding simulation for one AP is shown in Video S1. To see this figure in color, go online.

Video S1. Multi-scale Simulation

Time course of spatially detailed Ca2+ concentrations and the AP, related to Fig. 3.

The jSR concentration decreases rapidly upon the onset of release (Fig. 3). In case of sparks, this helps terminating release due to decreasing release current and consequently less coupling of RyRs by CICR. This mechanism is in agreement with earlier studies (25, 26, 44). In case of AP simulations, the jSR concentration continues to decrease on average till about the end of the membrane potential plateau (Fig. 3).

Fig. 3 also shows the cytosolic concentration [Ca2+]i caused by the release from all CRUs in the z-disk (15 μm × 15 μm × 2 μm) and the membrane potential, average [Ca2+]i, and average nSR concentration. The variable size of the volumes enclosed by the 2.6 μM isoconcentration surface illustrates the randomness and heterogeneity of release events. Refined numerical grids around CRUs guarantee the faithful simulation of Ca2+ and buffer diffusion between them (6). In that way, we can simulate the concentration dynamics from subdyadic to cellular length and timescales.

To facilitate the comparison of experimentally measured and simulated [Ca2+], we simulated a fluorescent buffer (see Fig. S4). This allowed us to emulate the approximation of [Ca2+]i as it would be measured by a single wavelength Fluo-4 experimental recording using an in vitro calibration approach as described in (45):where K is the dissociation constant of Fluo-4, F is the experimentally measured fluorescence intensity (the spatial average of b), Fmax is the measured fluorescence intensity in Ca2+-saturated dye (here, this is set as ), and Fmin is the measured fluorescence intensity in the absence of Ca2+ (here, set to zero).

Quantification of parameter values based on biomarkers

A population of simulated cells was generated as described in Materials and Methods. We identified valid parameter sets by filtering all simulation results for those providing biomarker values in the ranges stated in literature (see Tables 2 and 3). The biomarker resting membrane potential, maximum membrane potential, dome membrane potential, and [Na+] (see Table 2) are mainly determined and met by the Mahajan electrophysiology model we use and have essentially not been affected by the parameter variations considered here (36). Results of AP simulations were filtered by taking the biomarkers APD90, peak systolic [Ca2+], and mean diastolic [Ca2+] into account. The spark biomarkers used in the filtering of the spark simulations are the spark rate (i.e., events with at least two simultaneously open RyRs in the same cleft), the average FDHM, the mean of the Ca2+ peak value of sparks, and the ratio between quarks and sparks. Quarks are events in which exactly one RyR opens in a given cleft. The overlap of both filtering results led to the suggested parameter value set in Table 1. We performed 281 simulations for APs. Of those, 23 parameter sets passed the biomarker ranges stated above. Out of the 297 Ca2+ spark simulations, 20 simulations passed the ranges for the Ca2+ sparks. The suggested value (Table 1, fourth column) fulfills the requirement for all seven biomarkers.

Table 2

Biomarkers Ranges for Action Potentials

Biomarker	Range	Description
Max Vm	46 ± 4.5 mV	maximal value of AP peaks
Resting Vm	−77.4 ± 3.9 mV	resting value of the AP
Dome Vm	15.2 ± 10.1 mV	peak in the plateau phase
APD90	150–200 ms	APD at 90%
Systolic [Ca2+]	0.6–1.2 μM	peak systolic calcium
Diastolic [Ca2+]	0.1–0.25 μM	diastolic calcium
[Na+]i	10.5–11.5 mM	intracellular sodium

Table S1 also lists references.

Table 3

Biomarker Ranges for Ca2+ Sparks

Biomarker	Range	Description
FDHM	8.0–17.5 ms	FDHM
Spark rate	1–5 μm−1s−1	number of sparks per μm cell and second
Quark to spark ratio	0.2–1.1	number of quarks/number of sparks
Peak [Ca2+]	10.0–22.0 μM	average maximal [Ca2+]
Peak [Ca^2+]iexp	0.6–1.2 μM	during a spark

Table S1 also lists references.

Sensitivity analysis and response surfaces of AP biomarkers

Sensitivity analysis provides information on how changes of input parameters affect a particular biomarker value. We have chosen to vary the five parameters kplus, kclose, gRyR, gLCC, and VP,max in this analysis. All of them are related to Ca2+ as the focus of this study. The parameters setting the RyR open probability (kplus, kclose) and SERCA uptake (VP,max) were chosen because they are targets of drugs or post-translational modifications. We vary the RyR conductivity gRyR because the in vivo single channel current is not well known. The LCC conductivity gLCC turned out to be an important parameter in preliminary simulations.

We have chosen the suggested values of the parameters in Table 1 as reference for all sensitivities in this study. Fig. 4 shows the sensitivities for the AP biomarkers (APD90, peak systolic [Ca2+], and mean diastolic [Ca2+]) versus the varied parameters. The LCC permeability gLCC has the strongest positive impact on all three biomarkers (Fig. 4). The positive correlation of APD90 with gLCC we observe is in line with results by Britton et al. (46). AP durations (APDs) are positively influenced by the RyR opening rate kplus and negatively by their closing rate kclose, whereas for the systolic [Ca2+]i peak and diastolic [Ca2+]i, the opposite is true. The influence of the opening and closing rates of the RyRs on the APD mediates the effect these two parameters have on the [Ca2+]i values. An increase of kplus decreases [Ca2+]i because it prolongs the AP. This entails longer Ca2+ release, which, in the end, reduces SR [Ca2+] and the release and leak currents. We see the opposite effect when increasing kclose. It shortens the AP and increases SR [Ca2+] and hence also the release and leak current.

Figure 4

Sensitivity values for APD90, systolic , and diastolic . Black bars represent the SD for the sensitivity coefficients obtained from Bayesian inference and indicate the uncertainty in the estimates. To see this figure in color, go online.

Extraction of more detailed information is based on the response surfaces. They provide an approximation for the dependency of each biomarker on the varied parameters. We use them to draw contour plots for biomarker values in dependency on parameter values (Fig. 5). The values of the parameters not varied in these plots are listed in Table 1 (fourth column) and Tables S2–S9.

The strong influence of gLCC on the AP biomarkers motivates the focus mainly on this parameter and how its change might be compensated for by a change of another parameter. Fig. 5 depicts biomarker values computed in the gLCC-VP,max plane. The iso lines for the upper and lower parameter values of the literature ranges are color coded in green and orange, respectively. From Fig. 5, we can read off how coordinated parameter changes can maintain important biomarkers as, for instance, APD90. Interestingly, the contours of iso-APD90 and isosystolic [Ca2+] are similar, and a coordinated change of SERCA uptake and LCC current along them could maintain approximately both but would affect mean diastolic [Ca2+].

Sensitivity analysis and response surfaces of spark biomarkers

Fig. 6 shows the sensitivities of the spark biomarkers with respect to the varied parameters. The spark rate is mainly influenced by the opening probability of the RyR kplus and the RyR permeability gRyR. Surprisingly, the full duration at half maximum (FDHM) is only weakly affected by all five parameters. Peak calcium strongly responds to changes of the RyR permeability. The quark to spark ratio is strongly negatively affected by changes in kplus and gRyR and positively by kclose.

Figure 6

Sensitivity values for spark rate, mean peak [Ca2+], spark FDHM, and quark to spark ratio. Black bars represent the SD for the sensitivity coefficients obtained from Bayesian inference and indicate the uncertainty in the estimates. To see this figure in color, go online.

The response surfaces of FDHM, spark rate, and are depicted in Fig. 7. was calculated from Ca2+-bound dye buffer, as in experimental analyses, by using Eq. 3. Although this inferred [Ca2+] has a reasonable accuracy for the mean [Ca2+] during APs (see Fig. 3), our simulations suggest that it fails for the quantification of spark peak [Ca2+]i (compare Figs. S4 and S7). Whereas the true peak [Ca2+]i reaches values of more than 20 μM for single sparks, the inferred experimental concentration only reaches values slightly above 1 μM. A similar discrepancy occurs for the spark FDHM. Although the true underlying spark events appear to be short release events with a duration of 5–15 ms, the FDHM for the inferred is by factor of ∼2 longer.

Functional consequences of geometrical properties of RyR clusters

We start with comparing two different models for RyR placement in clusters. The first model assumes a regular arrangement with equidistant spacing (40 nm) of RyRs on a regular grid, whereas the second one assumes irregular clustering properties on the basis of the measurements of Jayasinghe et al. (17) as described before. We use the suggested values of the parameter set in Table 1 and performed 10 simulations for each placement model. Channel numbers for the individual CRUs are drawn from the same distribution (Eq. S12) for both groups. Differences between the individual simulations even within one placement model group arise from the randomness of channel numbers in CRUs and their placement. We find clear differences of biomarker values between the two placement models (Figs. 8 and 9). Hence, the dyadic substructure clearly affects cellular responses.

Figure 8

AP biomarker values APD90, systolic and diastolic , and time to peak from simulations with two different RyR placement models, which are explained in the text. The boxplots are standard box and whisker diagrams. To see this figure in color, go online.

Figure 9

Spark biomarker values FDHM, peak , spark rate, and quark to spark ratio from simulations with two different RyR placement models. The boxplots are standard box and whisker diagrams. To see this figure in color, go online.

Geometric effects within the placement model by Jayasinghe et al. are shown in Fig. 10. A variety of channel configurations has been generated by sampling from the placement model distributions and additionally varying the distribution parameters. Although channel configurations have been characterized successfully by the adjacency matrix (22, 26), we are looking here for a simpler approach. We characterized channel configurations by a variety of measures (average nearest- and four nearest-neighbor distance, area per channel determined by convex hull, and mean occupancy (see Fig. S5)) and found mean occupancy to show the strongest correlation with biomarker values. Mean occupancy is 1 if all RyRs are far apart and 0 if all RyRs are in the same spot (Fig. S5).

Figure 10

The relation between mean occupancy, as described in Fig. S5 and averaged over all CRUs, and the biomarkers APD90, systolic and diastolic , and time to peak. Each dot corresponds to one simulation. Black dots mark Jayasinghe placement, and red circles mark regular placement. The Jayasinghe placement with the measured parameters has a mean occupancy of ∼0.49. The slopes of a linear regression for the Jayasinghe placement are as follows: APD90 −38 ms, systolic −0.68 μM, diastolic −0.056 μM, and time to peak 13 ms. To see this figure in color, go online.

APD90, peak systolic , and diastolic decrease with increasing mean occupancy (Fig. 10). This concerted decrease reflects the correlation between these values found in the contour plots in Fig. 5, too. All three trends are in line with the general picture of decreased Ca2+ release due to increased mean occupancy. Hence, we find effects of dyadic substructure also within one placement concept.

The regular placement exhibits smaller APD90, larger peak systolic , and shorter time to peak (Fig. 10). We assume that this is caused by the relation between LCC placement and the overall cluster shape. LCC locations were chosen according to the same rules described above for both placement methods. The Jayasinghe placement produces elongated clusters, and the regular placement produces quadratic clusters. Hence, the average RyR distance to the closest LCC is smaller with the regular placement than with the Jayasinghe placement, which entails stronger LCC-RyR coupling. This stronger and earlier Ca2+ release causes faster Ca2+-dependent inhibition of LCCs and thus shorter APD.

The strength of the coupling of RyRs by Ca2+ diffusion decreases with increasing mean occupancy. The averages of spark biomarker values depend in the expected manner on occupancy as the slopes of the linear regressions show (Fig. 11). FDHM increases with increasing mean occupancy, reflecting the known phenomenon of slower termination of sparks with weaker spatial coupling of RyRs. Correspondingly, peak systolic and spark rate decreases with increasing mean occupancy. Large quark to spark ratios were found with weaker RyR coupling at large mean occupancy only. The spark biomarker values exhibit much stronger fluctuations than the AP simulations. The results with regular placement fit into the relation on mean occupancy.

Figure 11

The relation between mean occupancy, as described in Fig. S5 and averaged over all CRUs, and the biomarkers FDHM, peak , spark rate, and quark to spark ratio. Each dot corresponds to one simulation. Black dots mark Jayasinghe placement, and red circles mark regular placement. The Jayasinghe placement with the measured parameters has a mean occupancy of ∼0.49. The slopes of a linear regression for the Jayasinghe placement are as follows: FDHM 7.42 ms, peak −0.88 μM, spark rate −0.37 s−1, and quark to spark ratio −0.11. To see this figure in color, go online.

Only the specific realizations of channel numbers in the individual CRUs vary between the individual simulations with the regular placement. Simulations with the Jayasinghe placement are distinguished by both channel number realizations and specific placement. Hence, comparing the scatter of the biomarker value results with regular placement (red circles) with the Jayasinghe placement (black dots) in Figs. 10 and 11 provides an idea of how much of the variability is due to the randomness of channel numbers per CRU. The variability due to channel number randomness is comparable to the total variability for the spark biomarkers FDHM and rate (Fig. 11). Jayasinghe placement increases the quark to spark ratio variability with increasing mean occupancy because coupling between channels becomes weaker. Surprisingly, the regular placement has a larger variability of the peak than the Jayasinghe placement because it exhibits also very large values. The same comparison for the AP simulations suggests variability of peak systolic to result mainly from channel number variability. Variability of APD90, diastolic , and time to peak increase substantially because of the Jayasinghe placement.

Discussion

Parameterization of detailed multiscale models faces the problem of large computational costs required for simulations, which often prevents systematic parameter searches. These models require large efforts to reduce compute time and methods concluding parameter dependencies from a minimum of simulations. We reported our efforts to speed up simulations in previous studies (1, 6). Here, we present a parameterization of the model to experimental data. To that end, we identified ranges of altogether 11 biomarker values from literature, four of them were fulfilled by model output, mainly based on previous work by Mahajan et al. (36), and seven of them were affected by our spatially detailed approach and had to be met by our parameterization procedure. The multiscale set up of our model allowed for detailed channel placement inside dyadic clefts, according to measurements in Jayasinghe et al. (17). Accordingly, we focused here on parameters that are crucial for the Ca2+ handling in the microscopic domain of the dyadic cleft: the opening and closing rates of RyRs kplus and kclose, the RyR permeability gRyR, the LCC permeability gLCC, and the strength of SERCA uptake VP,max.

In the case of AP simulations, LCC permeability gLCC variation has the strongest impact on cell behavior, as illustrated by the sensitivities in Fig. 4. The response surfaces in Fig. 5 illustrate the strong correlation between APD and the systolic [Ca2+]i values. As an example of a control strategy suggested by response surfaces, we note that a coordinated change of gLCC and VP,max can approximately maintain APD and systolic [Ca2+]i while lowering diastolic [Ca2+]i, if the initial APD is in the lower range of the literature values. That would be a strategy to decrease the propensity for diastolic triggered events while maintaining contraction.

The spark simulations show that kplus, kclose, and gRyR affect the spark rate in the way we expected: kplus and gRyR positively, kclose negatively. Their effect on the quark to spark ratio can be comprehended by considering the probability that the first open RyR does open another one, thus turning a quark into a spark. That probability increases with kplus and gRyR and therefore reduces the quark to spark ratio and vice versa with kclose. The weak effect of all five parameters on the FDHM is surprising at a first glance; however, it agrees with the results by Cannell et al. (25) for the latency of induction decay.

The Jayasinghe placement affects the mean and SD of the AP biomarker values and mean and/or the SD of the spark biomarker values we have investigated. Hence, subdyadic structure matters. Figs. 8, 9, 10, and 11 report the effect of the Jayasinghe placement compared to the regular placement. These effects are due to cluster shape and channel density. The elongated shape of the clusters generated by the Jayasinghe placement weakens LCC-RyR coupling (compared to the regular placement), which increases time to peak and APD90 and affects indirectly diastolic and peak systolic . Mean occupancy is related to density and mean channel distance. In principle, the SD of channel distances could also have an effect on CRU dynamics by generating highly coupled subclusters at large SD. However, that is not supported by our simulations.

We found the effects of channel placement, which might be affected by specifics of our cleft model. We assume a fraction of 50% of total Ca2+ to be buffered by mobile buffers based on estimates from (26). This estimate assumes ATP (K ≈ 200 μM (26)) to be the dominating mobile buffer at Ca2+ concentrations occurring in the dyadic space. ATP levels may change to a degree affecting dyadic buffering in pathological states; however, this was not in the scope of this study (47, 48). We use an effective diffusion coefficient of 100 μm2 s−1 inside the cleft as well as a quasistatic approximation for the concentration profiles (see Supporting Materials and Methods; (1, 6, 32, 33)). Coupling of RyRs upon the opening of a channel with dynamic concentration profiles is initially weaker than with quasistatic profiles. Geometric effects are more important with weak spatial coupling. Hence, we assume that dynamic profiles would slightly amplify them. The timescales of the cytosolic concentration around the cleft space dominate the dynamics upon closing. They are captured by our model. The diffusion coefficient in the dyadic space has not been measured, and we informed our model on the basis of other modeling studies (25, 26, 49). Because we are using an effective value, it is rather at the upper end of the currently accepted range. In summary, we may rather underestimate than overestimate the geometrical effects.

In this study, we used multiscale modeling in two ways complementary to experimental studies. We observed local dynamics and cellular behavior simultaneously, and we modified subdyadic structures, which cannot be modified by experimental means. That provided insight into structure-function relations across multiple scales. Our results suggest that both the mean occupancy and the overall cluster shape affect APs and cytosolic Ca2+ transients.

Author Contributions

M.F., W.G., F.G.C., N.C., U.P., and S.L. designed the research. W.G., F.G.C., and W.N. performed the research. W.G. and F.G.C. analyzed the data. M.F., F.G.C., W.G., N.C., U.P., and S.L. wrote the manuscript.