Optogenetic actuation in ChR2-transduced fibroblasts alter excitation-contraction coupling and mechano-electric feedback in coupled cardiomyocytes: a computational modeling study

1.   Introduction

Cardiac fibroblasts (Fbs), which are essential for maintaining homeostasis of cardiac structure and function under normal physiological conditions, can also lead to tissue remodeling to interfere with normal conduction of excitation in pathological wound healing ^[1,2]. Electrophysiological remodeling effects of Fbs on cardiomyocytes (CMs) including conduction slowing, action potential duration (APD) extending, resting membrane potential (${V}_{\text{rest}}$) increasing and upstroke velocity decreasing ^[3,4,5]. Studies also verify that Fbs can increase cardiac tissue contractility via calcium (Ca²⁺) transients, a direct determinant of excitation-contraction coupling (ECC) ^[6,7].

The roles of Fbs on CMs mentioned above are observed by conventional electrical and biochemical methods. Contemporaneously, cardiac optogenetics has been emerging as an exciting area involving the delivery of opsins to enable optical modulation of cardiac electrical function ^[8]. With noteworthy advantages over conventional methods, such as selective electrophysiological modulation and high-resolution spatiotemporal control, this approach has been widely applied in exploring whether cardiac Fbs couple electrically with working myocytes ^[9,10,11]. In the process of optogenetic stimulations in CMs and nonmyocytes, the light-sensitive ion channel, Channelrhodopsin-2 (ChR2), has emerged as an important excitatory optogenetic tools in cardiac electrophysiology ^[12,13,14].

The ECC, which is essential to maintain calcium homeostasis and cardiac function, has been well studied so far ^[15]. While the inverse process, mechano-electric feedback (MEF), has scarce data from human cardiac experiments and need further research. In MEF, mechanical effects can modulate cardiac electrical activity in complex ways ^[16]. The basic phenomena of MEF include the length dependence of isometric contractions and the load dependence of isotonic afterloaded contractions, and reveal as responses of calcium transients and action potentials (APs) to mechanical interventions ^[17,18].

Most of the studies mentioned above are completed by experiments, which actually have made great strides towards demonstrating the cardiac mechanisms and the feasibility of clinical applications. However, computational modeling, as an alternative option, can help accelerate experimental progress by offering a proving ground to assess whether these proposed approaches could be scaled up and used in human hearts. The aim of this study is to assess the influences of optogenetic actuation on ECC and MEF in myocardium, which are, to our best knowledge, barely explored. To achieve this goal, we present a preliminary combined model incorporating optogenetic capabilities in myocardial ECC and MEF. We model photokinetics of the light-sensitive ChR2 channel and represent ChR2 delivery via Fbs, incorporate it into the CM-Fb coupled model, and evaluate ECC and MEF performances in the coupled cells at various light durations (LDs) and light irradiances (LIs).

2.   Materials and methods

Our combined model of the human CM coupled with ChR2-transduced Fbs is formed of three main modules: a module of ECC in the CM, a module of optical actuation in Fbs, and a module of CM-Fb coupling. A complete description of this model can be found in the Supporting Material.

2.1.   The module of ECC in the CM

This module contains two parts: an electrophysiological and a mechanical part coupled to each other.

The part for cellular electrophysiology is based on the 'Maleckar-Trayanova' (MT) ionic model, which is described in detail elsewhere ^[19]. This model describes ionic currents across the membrane and correctly replicates APD restitution of the adult human atrial myocyte. It has been widely used in electrophysiological studies and is used as the electrical component in electromechanical models ^[20,21,22].

The part of mechanical activity is based on the mechanical block of an updated version of the 'Ekaterinburg-Oxford' (EO) model, which is reported in detail elsewhere ^[23]. It consists of a contractile element, three elastic and two viscous elements. In particular, this model pays special attention to the processes of association and dissociation of calcium with troponin C (CaTnC) complexes along actin filaments, describes the length-dependent probability of cross-bridge (Xb) binding in detail, and thus reproduces most of ECC and MEF effects found in the heart muscle ^[17,24].

The above two parts are coupled to form the module of ECC in the CM via Ca²⁺ handling. The MT model uses fluid compartment formulations to describe the buffering action of calsequestrin, troponin and calmodulin buffers. The EO model comprises equations describing the kinetics of calcium complex formation and links them to Xb kinetics and sarcomere length. We replace the description of intracellular Ca²⁺ buffering in the MT model with that in the EO model, modify the formulation of intracellular Ca²⁺ concentration in the modified MT model, and then integrate the two models to compose the module of ECC. Here we provide only the equations of kinetics of CaTnC complexes and intracellular Ca²⁺ buffering that are used and modified in our ECC model.

The following equation describes time-dependent changes in the concentration of CaTnC complexes ($\left[\text{CaTnC}\right]$) ^[23]:

where ${a}_{\text{on}}$ is the rate constant for CaTnC association; ${\text{TnC}}_{\text{tot}}$ is the total concentration of troponin C (TnC); ${\left[{\text{Ca}}^{\text{2+}}\right]}_{i}$ is the intracellular Ca²⁺ concentration; ${a}_{\text{off}}$ is the maximum rate constant for CaTnC dissociation; ${k}_{\text{A}}$ is the cooperativity parameter; $\varPi \left({N}_{\mathrm{A}}\right)$ is the cooperative dependence defining cooperativity of the contractile proteins; ${N}_{\mathrm{A}}$ is the average function of the attached Xbs per one CaTnC complexes; ${\varPi }_{\text{min}}$ is the minimum value of the cooperative dependence; $N$ is the fraction of force generating attached Xbs; ${L}_{\text{oz}}$ is the instantaneous length of thick and thin filament overlap zone.

The following equations describe time-dependent changes in Ca²⁺ bound concentration and in free intracellular Ca²⁺ concentration ^[25,26]:

where ${B}_{1}$ and ${B}_{2}$ are the concentrations of Ca²⁺ bound with "fast" and ''slow'' Ca²⁺ binding ligands; ${b}_{\text{1on}}$ and ${b}_{\text{2}\text{on}}$ are the rate constants for Ca²⁺ binding; ${b}_{\text{1off}}$ and ${b}_{\text{2}\text{off}}$ are the max rate constants for Ca²⁺ unbinding; ${B}_{\text{1tot}}$ and ${B}_{\text{2}\text{tot}}$ are the total concentrations of Ca²⁺ binding ligands; ${I}_{\text{di}}$ is the Ca²⁺ diffusion current from the diffusion-restricted subsarcolemmal space to the cytosol; ${I}_{\text{B, Ca}}$ is the background Ca²⁺ current; ${I}_{\text{CaP}}$ is the sarcolemmal Ca²⁺ pump current; ${I}_{\text{NaCa}}$ is the Na⁺-Ca²⁺ exchange current; ${I}_{\text{up}}$ is the sarcoplasmic reticulum Ca²⁺ uptake current; ${I}_{\text{rel}}$ is the sarcoplasmic reticulum Ca²⁺ release current; ${\text{Vol}}_{\text{i}}$ is the total cytosolic volume; $\text{F}$ is the Faraday's constant.

2.2.   The module of optical actuation in Fbs

This module contains an electrophysiological model of Fbs transduced with ChR2 currents. In other words, the modular form of the ChR2 current model is direct plug-in into the electrophysiological model of Fbs. Current literature shows that optogenetic stimulation do not lead to changes in Fb force and strain generation ^[11], therefore, we don't add the mechanical model of Fbs into this module.

The Fb model used in this study is the human atrial Fb model developed by Maleckar et al., which is based on the general formulation of MacCannell et al. and modified to represent atrial Fb properties, and the Williams model for ChR2 current is added to it ^[12,27,28]. Here we provide only the current output for the ChR2-transduced Fb model (${I}_{\text{F}\text{b}}$).

where ${I}_{\text{f}}$ represents the transmembrane current across Fb; ${I}_{\text{ChR2}}$ represents the ChR2 current.

2.3.   The module of CM-Fb coupling

In this module, a selected number of Fbs are coupled to a CM by assigning an intercellular Fb-to-CM conductance (${G}_{\text{gap}}$). The following equations describe time-dependent changes in the transmembrane potentials of CM and Fb ^[27]:

where ${V}_{\text{Fb}i}$ represents the transmembrane potential across the $i$th coupled Fb; ${I}_{\text{F}\text{b}i}$ represents the current output for the $i$th coupled Fb; ${V}_{\text{M}}$ represents the CM membrane potential; ${\mathrm{C}}_{\text{M}}$ and ${\mathrm{C}}_{\text{Fb}}$ represent the CM and Fb membrane capacitances, respectively; ${I}_{\text{M}}$ represents the CM transmembrane current as defined by the MT mathematical model ^[19]; $n$ is the total number of Fbs coupled to the CM.

2.4.   Simulation protocol

The prime objective of this manuscript contains two parts. One is to analyze possible roles of optogenetic actuation on the ECC of CM and the electrophysiology of Fb at various LDs and LIs. The other is to compare the MEF effects of CM at specified values of LD and LIs. our research is focused on these effects caused by light rather than CM-Fb coupling, hence, parameters related to CM-Fb coupling, such as $n$ and ${G}_{\text{gap}}$, are fixed in our simulations ($n = \text{2}$, ${G}_{\text{gap}} = 2.0 \ \ \text{n}\text{S}$).

The combined model is optically paced at a rate of 1 Hz. To ensure the system reached a steady-state, simulation results of CM and Fb electrophysiology are recorded after 100 optical stimuli, and simulated responses to mechanical interventions (such as changes in myocyte length and afterload) are evaluated during the first twitch in which the intervention is applied according to experimental protocols ^[29,30].

The coupled model is implemented in MATLAB with a robust protocol for optical stimulation. Integration for the cell models is done using a built-in integration algorithm with a variable time step and suited for stiff systems of ordinary differential equations (ode15s) at absolute and relative error tolerances of ${\text{10}}^{\text{-10}}$.

3.   Results

3.1.   Effects of optogenetic actuation on CM ECC and Fb electrophysiology at various LDs and LIs

Figure 1 shows the ECC in CM and the electrophysiology in Fb triggered by a light pulse with a LD of 50 ms and a LI of 5 mW/mm².

The top panel shows the optically-triggered ${I}_{\text{ChR2}}$ in Fbs and membrane potentials in the CM and Fbs. Upon a typical supra-threshold light pulse, ChR2 channels are opened and an inward ${I}_{\text{ChR2}}$ is generated, which causes an AP like waveform in Fbs. Changes in the membrane potential in Fbs further transmit to the CM via gap junction, elicit CM depolarization and finally initiate an AP in the CM.

During the excitement of CM, the kinetics of CaTnC complexes regulates the CM's mechanical activities. The middle and bottom panels illustrate the dynamic changes in force ($F$), muscle and sarcomere length ($L$ and $SL$) in isometric and isotonic contractions. In the simulation of isometric twitch, the constant muscle length, increased isometric force and shortened internal sarcomere are observed, which match the experimental data obtained in isolated human trabeculae ^[31]. During the isotonic afterloaded contraction, the lengths of muscle and sarcomere decrease, which also fits well into the range of experimental data ^[31,32].

Figures 2 and 3 show changes in ${V}_{\text{rest}}$ and the amplitude of the membrane potential (${V}_{\text{max}}$) of the CM and Fbs as a function of LD and LI. The APD at 80% repolarization (${APD}_{\text{80}}$) of the CM is also recorded. In our simulations, the inward current ${I}_{\text{ChR2}}$ in Fbs causes an AP like waveform in ${V}_{\text{Fb}}$. However, "APD" is not applicable to Fbs because they are rather passive and cannot initiate regenerative APs as the CM. Therefore, we rephrase it as "membrane potential duration at 80% restitution (${MPD}_{\text{80}}$)" for Fbs.

In Figure 2, LI is fixed as 5 mW/mm² and LD varies from 30 to 210 ms with an increment of 10 ms. In general, the ${V}_{\text{rest}}$ of both cells are lowered and then raised as LD increases. The maximum and minimum of ${V}_{\text{rest}}$ in the CM and Fbs are −67.7 mV vs −68.5 mV and −65.5 mV vs −66.3 mV, respectively. As LD increases from 30 ms to 50 ms, the ${V}_{\text{max}}$ in the CM and Fbs are raised from 13.3 to 14.9 mV and from 2.7 to 4.9 mV, respectively, and keep stable before LD is 170 ms. When LD continues to grow, small fluctuations appear in ${V}_{\text{max}}$. ${APD}_{\text{80}}$ and ${MPD}_{\text{80}}$ drop slightly at first and increase steadily when LD is greater than 50 ms.

In Figure 3, LD is fixed as 50 ms and LI varies from 0.3 to 10 mW/mm² with an increment of 0.1 mW/mm². The ${V}_{\text{rest}}$ of both cells decrease gradually as LI varies from 0.3 to 0.9 mW/mm², and keep at stable levels (~−68.1 mV for the CM and ~−65.9 mV for Fbs) when LI continues to grow. The trends of ${V}_{\text{max}}$ are almost opposite of ${V}_{\text{rest}}$ as LI is enhanced from 0.5 to 10 mW/mm². They increase gradually and then keep at ~14 mV for the CM and ~5.5 mV for Fbs when LI is greater than 1.1 mW/mm². When LI is 0.3 or 0.4 mW/mm², ${V}_{\text{max}}$ of the CM are −49.3 or −46.4 mV, which leads to failed excitation in the CM. Hence, there are no ${APD}_{\text{80}}$ and ${MPD}_{\text{80}}$ at these two LIs. ${APD}_{\text{80}}$ and ${MPD}_{\text{80}}$ are recorded when LI is between 0.5 and 10 mW/mm². The trends of ${APD}_{\text{80}}$ and ${MPD}_{\text{80}}$ are similar to ${V}_{\text{rest}}$. Drop at the beginning and keep steady at ~99 ms for the CM and ~91 ms for Fbs.

In Figure 4, a broader parameter space is explored, where LD and LI are varied from 30 to 210 ms and 0.3 to 10 mW/mm², respectively. When both LD and LI are very low (the bottom left corner of Figure 4A, 4B, 4D and 4E), The ${V}_{\text{rest}}$ and ${V}_{\text{max}}$ of both cells are at their highest and lowest levels, respectively. Outside these areas, changing LD or LI does not bring significant changes in ${V}_{\text{rest}}$ and ${V}_{\text{max}}$ of both cells. For each level of LI, ${V}_{\text{rest}}$ in both cells are lowered and then raised as LD prolongs. For each level of LD, it remains roughly unchanged as LI enhances. When LI increases from 0.3 to 2 mW/mm², the increments of ${V}_{\text{max}}$ of both cells are more obvious than the increments induced by the growth of LDs. The model shows that Fbs have no ${MPD}_{\text{80}}$ and the CM remains inexcitable at low LDs and LIs (Figure 4C and 4F, the black region), since either there is too little ChR2 current produced, or this current is unable to change Fbs membrane potentials apparently and to depolarize CMs, respectively. For the gray region in Figure 4C and 4F, ${MPD}_{\text{80}}$ and ${APD}_{\text{80}}$ are not calculated because spontaneous beatings are produced and exceed the pacing rate. The region is located in a LD between 190 and 210 ms and a LI between 1.4 and 3.5 mW/mm². Except the inexcitable and spontaneous beating regions, an increased LD prolongs ${MPD}_{\text{80}}$ and ${APD}_{\text{80}}$ with fixed LIs, while an enhanced LI does not invoke significant changes in them with fixed LDs.

Figure 5 is the simulation of isometric twitches (top panel) and isotonic afterloaded twitches (bottom panel) of CM obtained at various LDs (increased from 30 to 210 ms with an increment of 45 ms). LI is fixed to 5 mW/mm². The afterload applied in isotonic twitches is 3 mN. When LD is prolonged, the virtual muscle responds with an increase in peak force and twitch duration. The time to peak force is also increased from 268 to 275 ms (Figure 5A). The muscle length during isometric contractions keeps constant (Figure 5B). Internal sarcomere shortening is observed. Longer LD increases the duration of sarcomere deformation (Figure 5C). For isotonic twitches, as the LD prolongs, the time to plateau force is decreased slightly from 103 to 95 ms and the duration of plateau force is extended from 328 to 356 ms (Figure 5D), the end-systolic lengths of both muscle and sarcomeres are decreased while the deformation durations are increased. Meanwhile, the time to valley $L$ and $SL$ is increased from 305 to 317 ms (Figure 5E and 5F).

Figure 6 illustrates the effects of LI on isometric twitches (top panel) and isotonic afterloaded twitches (bottom panel). LI increases from 0.3 to 10 mW/mm² with an increment of 2.4 mW/mm². LD is fixed to 50 ms. The afterload in isotonic twitches is also 3 mN. Compared to the influence of LD on mechanical properties, the impact of LI is much less. When LI is 0.3 mW/mm², the ChR2 current is too little to depolarize the CM, and therefore, no twitch takes place. The force and lengths don't change. At other four LIs, changes in each parameter in both isometric and isotonic contractions are substantially the same (Figure 6A to 6F).

3.2.   Effects of optogenetic current on CM MEF at specified LDs and LIs

We also evaluate MEF manifestations in the CM when LD and LI are set to be 30 or 210 ms, and 1 or 10 mW/mm², respectively.

In Figure 7, we simulate isometric twitches with various diastolic lengths of the virtual preparation at a LD of 30 ms and 210 ms. LI is fixed as 5 mW/mm². The initial length (${L}_{\text{init}}$) of the sample is decreased from 95% to 75%${L}_{\text{max}}$ (${L}_{\text{max}}$ corresponds to a sarcomere length equal to 2.3 m). Both mechanical and electrical readouts are presented.

For both LDs, our simulation shows: (1) the virtual muscle responds with a decrease in the peak isometric $F$ when stretched less (Figure 7A and 7G); (2) length-dependent changes in the active force are related to a decrease in initial $SL$ and an increase in end-systolic sarcomere shortening (Figure 7B: 0.234 m per sarcomere at 95%${L}_{\text{max}}$ versus 0.245 m at 75%${L}_{\text{max}}$; Figure 7H: 0.238 m per sarcomere at 95%${L}_{\text{max}}$ versus 0.249 m at 75%${L}_{\text{max}}$); (3) CaTnC binding has a lower peak and decreases quickly at smaller $SL$ (Figure 7D and 7J). No obvious effect of stretch on ${V}_{\text{M}}$, ${\left[{\text{Ca}}^{\text{2+}}\right]}_{i}$ and Na⁺-Ca²⁺ exchange current (${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$) is found (Figure 7C, 7E, 7F, 7I, 7K and 7L).

Compare the two panels of Figure 7, we can see the effects of LD on CM MEF (from Figure 7A vs 7G to 7F vs 7L). For each ${L}_{\text{init}}$, as LD increases, (1) durations of twitch, sarcomere deformation, membrane potential and $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$ are extended, and ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ returns to the reverse mode again; and (2) peaks of $F$, ${V}_{\text{M}}$ and ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ increase slightly (~0.56 mN↑, ~3.78 mV↑ and ~0.73 pA/pF↑, respectively) while peaks of $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$ and ${\left[{\text{Ca}}^{\text{2+}}\right]}_{i}$ and minimums of $SL$ decrease slightly (~0.07 μM↓, ~0.004 μM↓ and 0.004 μm↓, respectively).

Figure 8 shows isotonic twitches with various afterloads at a LD of 30 ms and 210 ms. LI is fixed as 5 mW/mm². The afterload is decreased from 90%${F}_{\text{isom}}$ to 10%${F}_{\text{isom}}$ (${F}_{\text{isom}}$ is a peak active isometric force at ${L}_{\text{init}}$).

For both LDs, our simulation shows: (1) the duration of afterloaded twitches shortens with decreasing afterload (Figure 8A and 8H); (2) $L$ and $SL$ decrease, while the velocities of shortening and lengthening increase as the afterload is reduced (Figure 8B, 8C, 8I and 8J); (3) CaTnC complexes dissociate becomes fast with a decrease in afterload (Figure 8E and 8L), which results in a prolongation of Ca²⁺ transients (Figure 8F and 8M), and accelerates ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ to flow inward (Figure 8G and 8N) and further prolongs APD (Figure 8D and 8K).

Compare the two panels of Figure 8 (from Figure 8A vs 8H to 8G vs 8N), we can see the effects of LD on CM MEF are similar to that in Figure 7. For each ${F}_{\text{afterload}}$, as LD increases, (1) durations of twitch, muscle and sarcomere deformation, membrane potential, Ca²⁺ transient and ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ are extended; (2) peaks of ${V}_{\text{M}}$ and ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ increase slightly (~3.78 mV↑ and ~0.73 pA/pF↑, respectively) while peaks of $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$ and ${\left[{\text{Ca}}^{\text{2+}}\right]}_{i}$ decrease slightly (~0.06 μM↓ and ~0.005 μM↓, respectively); and (3) end-systolic lengths for both muscle and sarcomeres decrease.

Figure 9 shows isometric twitches with various initial sample lengths at a LI of 1 mW/mm² and 10 mW/mm². The LD is fixed as 50 ms (since the CM cannot be excited by a LI of 0.3 mW/mm² when LD is 50 ms, here we choose a larger LI of 1 mW/mm²). ${L}_{\text{init}}$ of the sample is decreased from 95% to 75%${L}_{\text{max}}$.

For both LIs, the MEF effects at each ${L}_{\text{init}}$ showed in Figure 9 (9A to 9F, and 9G to 9L) are similar to that in Figure 7. While the effects of LI on CM MEF (from Figure 9A vs 9G to 9F vs 9L) are not the same as that of LD (Figure 7). For each ${L}_{\text{init}}$, as LI increases, (1) durations of twitch, sarcomere deformation, membrane potential, Ca²⁺ transient and ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ are unchanged; (2) peaks of ${V}_{\text{M}}$ and ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ increase slightly (~6.76 mV↑ and ~1.33 pA/pF↑, respectively) while peaks of $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$ and ${\left[{\text{Ca}}^{\text{2+}}\right]}_{i}$ decrease slightly (~0.28 μM↓ and ~0.01 μM↓, respectively); and (3) peaks of $F$ and minimums of $SL$ are barely changed (~0.1 mN↓ and ~0.00007 μm↓, respectively).

Figure 10 shows isotonic twitches with various afterloads at a LI of 1 mW/mm² and 10 mW/mm². The LD is fixed as 50 ms. The afterload is decreased from 90%${F}_{\text{isom}}$ to 10%${F}_{\text{isom}}$.

For both LIs, the MEF effects at each ${F}_{\text{afterload}}$ showed in Figure 10 (10A to 10G, and 10H to 10N) are also similar to that in Figure 8. While the effects of LI on CM MEF (from Figure 10A vs 10H to 10G vs 10N) are not the same as that of LD (Figure 8). For each ${F}_{\text{afterload}}$, as LI increases, (1) durations of twitch, muscle and sarcomere deformation, membrane potential, Ca²⁺ transient and ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ are unchanged; (2) peaks of ${V}_{\text{M}}$ and ${I}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}$ increase slightly (~6.76 mV↑ and ~1.33 pA/pF↑, respectively) while peaks of $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$ and ${\left[{Ca}^{2+}\right]}_{i}$ decrease slightly (~0.28 μM↓ and ~0.01 μM↓, respectively); and (3) minimums of $L$ and $SL$ are barely changed (~0.1%↓ and ~0.004 μm↓, respectively).

4.   Discussion

From the modeling standpoint, the approaches used to describe light sensitivity at the myocyte level are generally classified into two categories: (1) direct expression of opsins in CMs and (2) electrical coupling of exogenous opsin-expressing cells to CMs ^[33]. Our modeling applies the second approach. Experimental works have confirmed that opsin-transfected human cells (e.g., embryonic kidney cells, NIH-3T3 Fb cells, cardiac Fbs and stem cell derived CMs) provide light sensitivity when they are co-cultured with CMs or injected into the beating muscle ^{[13,34,35,36]}. In this study, we present a combined model for incorporating optogenetic capabilities in detailed simulations of myocardial electrophysiology and mechanics. We model photo kinetics of the light-sensitive ChR2 channel and represent ChR2 delivery via Fbs, and then we incorporate it into the model of CM coupled with Fbs, and evaluate ECC and MEF manifestations in the coupled model at various LDs and LIs. We illustrate the optical capability from three aspects: first, we show that the increased LD and LI influence potential characteristics of the CM and Fbs; second, we show that the increased LD and LI influence mechanical activities in the CM; Finally, we show MEF performances in the CM at specified LDs and LIs. Here, we discuss results presented in our study.

Experimental observations have demonstrated that the minimal coupling conductance between ChR2-rich somatic cells and CMs required for cell delivery mediated optogenetics is approximately 2 nS ^[13]. In our simulations, ${G}_{gap}$ is set to be this value. The characteristics observed in the top panel of Figure 1 qualitatively resembles the experimentally observed response to optical stimulation in tandem preparations of ChR2-rich HEK cells and myocytes ^[13]. Meanwhile, we compare a direct electrical stimulation in a CM and an optogenetic stimulation in Fbs to excite the coupled CM indirectly (Figure S1 in the Supplementary Material). It shows that both short rectangular electrical pulse and long optical pulse can excite the CM, which is accordant with the findings by Williams et al. ^[37]. For CM-Fb coupling, the traditional electrical method can only excite CMs but is useless to Fbs, while light can induce steady inward current in ChR2-tranduced Fbs to regulate their membrane potential and further depolarize the coupled CM.

We also compare the intracellular Ca²⁺ concentration and AP waveform generated by the combined model with the original MT model at a rate of 1 Hz (Figure S2). The two models are in the same electrical stimulating environment and are not coupled with Fbs. The only difference is the description of Ca²⁺ handling. In our formulation, peak ${\left[{\text{Ca}}^{\text{2+}}\right]}_{i}$ is lower and the decrease of ${\left[{\text{Ca}}^{\text{2+}}\right]}_{i}$ is faster than in the original model. The AP waveform in the combined model has a slight decreased rate of late repolarization compared with the original model.

There has been an awful lot of good work done on light sensitization. Abilez et al. assess electrical, biochemical, and mechanical signals in ChR2-expressing CMs at three light stimulation frequencies. In their records, the evoked signals during light stimulation are markedly different from pre- and poststimulation signals at all frequencies ^[14]. Williams et al. demonstrate that ${I}_{\text{ChR2}}$ allows for optical control of APD in atrial and ventricular myocytes ^[37]. Karathanos et al. illustrate that abnormal APD shortening could be almost completely eliminated by illuminating ChR2-rich cells with appropriately shaped optical pulses ^[38]. Kostecki et al. show that Fb-specific optogenetic induced inward currents decrease conduction velocity, shorten APD and lead to spontaneous beatings in co-cultured CMs as LI increases ^[11]. In this study, we explore the roles of both LD and LI on the excitation and contraction of cardiac cells. Since many studies have demonstrated that CM-Fb coupling contributes to arrhythmia formation, we fix parameters related to this coupling (i.e., ${G}_{gap}$ and the number of coupled Fbs) to only focus on the electrophysiological and mechanical manifestations caused by light. In our simulations (Figures 2 to 4), a low LD and LI could not invoke the CM (for example, LD ≤ 40 ms and LI ≤ 1.5 mW/mm²). Further increase of LD and LI has a minor effect on ${V}_{\text{rest}}$ and ${V}_{\text{max}}$ of both cells, but a major effect on APD and MPD. On the one hand, spontaneous beatings are produced when LD is large and LI is relatively low (the gray region in the area of 190 ms ≤ LD ≤ 210 ms and 1.4 mW/mm² ≤ LI ≤ 3.5 mW/mm² in Figure 4C and 4F). Kostecki et al. have also verified spontaneous beatings in their CM-Fb coupled experiment and modeling, while they use a longer LD and a lower LI (LD = 2000 ms and LI between 0.057 and 1.71 mW/mm²) ^[11]. Compare their results to ours, it can be predicted that long low-intensity optical pulses could put a coupled tissue more at risk of spontaneous beatings and rhythm disorders. On the other hand, extreme low LDs and LIs cannot evoke the coupled cells (black regions in Figure 4C and 4F). Outside of these areas, the durations of membrane potential of both cells are extended as the LD prolongs, and remain relatively stable after initial short-lived fluctuations as LI enhances.

Since the initial findings of Bruegmann and Arrenberg et al. that optical stimulation could modulate pacemaker activities in transgenic mice and zebrafish heart, the research area about cardiac optogenetics is flourishing ^[39,40]. Until now, this technology has been transforming in the realm of cardiac electrophysiology research. In our study, not only the optical stimulation affects the membrane potential of CM (Figure 2 to 4), but also influences mechanical twitches by ECC (Figures 5 and 6) and in turn regulates electrophysiological properties of CM by MEF (Figures 7 to 10). Here we discuss the influences from three aspects.

The first is roles of LD and LI on CM's APD and twitch in the process of ECC. Figure S3 shows electrophysiological and mechanical parameters (in isometric contraction) related to ECC at a LD of 30, 120 and 210 ms. LI is fixed to 5 mW/mm². From the figure, we can see that the longest LD corresponds to the longest APD. Currents which may influence APD are surveyed (Figure S3F to S3J). The early phase of ${I}_{\text{ChR2}}$ is a short downward pulse. As LD increases, the inflow of ${I}_{\text{ChR2}}$ becomes strong again. This inward ${I}_{\text{ChR2}}$ during the repolarization phase is the main determining factor affecting APD prolongation. The influence of ${I}_{\text{ChR2}}$ on ${V}_{\text{M}}$ is also reflected by the curves of ${I}_{gap}$. After the early transient outward current, ${\text{I}}_{gap}$ changes from an outward current to an inward current when LD is increased, which tends to lengthen APD (the double nature of ${I}_{gap}$ on APD is well recognized by Xie et al. ^[41]). For the transient outward K⁺ current (${I}_{\text{t}}$) and ultra-rapid delayed rectifier K⁺ current (${I}_{\text{Kur}}$), their upstrokes are mostly unchanged as LD increases. The action potential at a longer LD achieves a higher level at the beginning of the plateau, thus the voltage-dependent relative changes in the inward rectifier K⁺ current (${I}_{\text{K1}}$) are enhanced, which also tends to prolong APD. Meanwhile, as LD increases, the dissociation of CaTnC complexes decreases, which further extends the duration of Xb attachment-detachment ($N$, a key parameter connecting CaTnC to $F$ and $SL$ in ECC, see Model equations in Supplementary Materials for more detail), finally increases peak force and durations of twitch and sarcomere deformation. Figure S4 shows parameters at a LI of 0.3, 5 and 10 mW/mm². LD is fixed to 50 ms. From this figure we can see that APD and the duration of twitch are barely changed when LI is enhanced. Due to the inward rectification kinetic feature of ChR2, AP response to optogenetic approach is slower than that to electrical stimulation, and the AP upstroke starts really later after significant delay (after the duration of illumination, as shown in Figure S1). As a result, the magnitude of the inward current becomes smaller at a higher voltage. It does not affect the membrane potential of Fb when ${V}_{\text{Fb}}$ is greater than 0 mV and its effect on the membrane potential of CM is limited when ${V}_{\text{M}}$ is above 0 mV. Since the AP upstroke above 0 mV is not affected by ${I}_{\text{ChR2}}$, kinetics of the voltage-dependent ${I}_{\text{t}}$ are not altered by this current. Figures S3 and S4 illustrate that the inward ${I}_{\text{ChR2}}$ and its inward rectification feature have dominant roles on CM's APD and twitch, and LD has dominant role than LI.

The second is roles of LD and LI on CM's twitch and APD in the reverse process, MEF. Figure S5 illustrates electrophysiological and mechanical parameters (in isometric and isotonic contractions) related to MEF at a LD of 30, 120 and 210 ms. LI is fixed to 5 mW/mm². In both twitches, LD prolongation increases peak force and twitch duration. This enhanced contraction intensifies internal sarcomere shortening and slows its relaxation, which delays the dissociation of CaTnC complexes. For the Na⁺-Ca²⁺ exchange, it acts in the reverse mode at first, and quickly switches to the forward mode after the AP upstroke. A major difference between the three Na⁺-Ca²⁺ exchange currents is that the current will return to the reverse mode again when LD is increased (see the curves of ${I}_{\text{NaCa}}$ that above 0 pA/pF between 146 and 227 ms in Figure S5D and between 171 and 220 ms in Figure S5I when LD is 210 ms). As LD increases, both CaTnC and Na⁺-Ca²⁺ exchange begin to contribute to against APD prolongation. However, relative changes in ${I}_{\text{ChR2}}$ tended to lengthen AP are so pronounced that oppositely directed relative changes in $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$ and ${I}_{\text{NaCa}}$ can no longer seriously influence the slowing of AP prolongation. Figure S6 illustrates parameters at a LI of 0.3, 5 and 10 mW/mm². The LD is fixed to 50 ms. This figure also shows that the duration of AP and twitch are barely changed when the CM is excited and LI continues to grow. As discussed above, the rectification feature of ${I}_{\text{ChR2}}$ limits its effect on myocyte membrane potential when ${V}_{\text{M}}$ is above 0 mV. Enhancing LI does not change this process and therefore cannot bring changes in APD. Meanwhile, $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$ and ${I}_{\text{NaCa}}$ are barely changed at the fixed LD and do not shorten AP.

The third is roles of ${L}_{\text{init}}$ and ${F}_{\text{afterload}}$ on CM's twitch and APD in MEF when LD and LI are fixed (each panel in Figure 7 to 10). Our results are in accord with previous modeling studies ^[17,42,43] and fit well into the range of experimental data obtained from animals ^[29,44]. Due to the lack of experimental data for the human CM, these results could not be directly confirmed or rejected by experiments so far. Balakina-Vikulova et al. have discussed the possible mechanisms in the above process in detail. They rule out the possibility of the involvement of stretch-activated channels (SACs) in these load- and length-dependent events, and indicate that the key factor which facilitate the above process is the CaTnC dissociation ^[17].

In MEF itself, it is revealed as length dependence of isometric twitches and load dependence of afterloaded isotonic twitches, and as respective reactions of APs to mechanical involvement. For each fixed LD or LI, MEF is found to have visible effects on the duration of $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$, ${I}_{\text{NaCa}}$ and AP in response to smaller afterloads (Figures 8 and 10). While if we compare the performances of these parameters related to MEF in conditions of an increased LD or LI and a constant ${L}_{\text{init}}$ or ${F}_{\text{afterload}}$, we can find that the change in ${I}_{\text{ChR2}}$ working to AP prolongation is more pronounced than changes in $\left[\mathrm{C}\mathrm{a}\mathrm{T}\mathrm{n}\mathrm{C}\right]$ and Na⁺-Ca²⁺ exchange contributing to AP shortening (Figures S5 and S6). Hence, the inward ${I}_{\text{ChR2}}$ and its inward rectification feature have dominant roles on CM's APD and twitch in both ECC and MEF.

Three limitations in this study should be mentioned. First, the MEF effects explored in our study are limited to fast responses of the CM to mechanical impacts (isometric and isotonic afterloaded ones). The slow responses of force and calcium, however, have not been evaluated. For example, the SACs, which are confirmed to present in CMs and regarded as a main factor underlying the slow response, may mediate Ca²⁺ entry during this process ^[45,46]. Changes in protein's expression and structure during the long-term MEF are also not dealt with in this study. Due to a large variation in reports of their effects on AP form and duration ^[47], MEF modeling including SACs should make a detailed description of precondition (such as the magnitude of the reversal potential and conductance of SACs, and the moment at which the mechanical impact is applied) to demonstrate the contribution of SACs to the development of AP. Second, only ChR2 is considered in this model. Other opsins (such as the light-sensitive chloride pump halorhodopsin and the proton pump Archaerhodopsin-T), proved to silence excitation ^[34,48], are not be characterized mathematically here. Third, the breadth of this computational study needs to be extended. Our work focused on the scale of local cell–cell interactions. Other scales, such as scales of cell–matrix interactions, tissue and organ level conduction properties, are not included in this preliminary study. In these scales, not only the effects of LD and LI, but also the influences of realistic spatial patterns of opsin-expressing cells and of the attenuation of light energy on cardiac electrophysiological and mechanical activities should be considered into optogenetics modeling.

5.   Conclusions

Here we use computational modeling to show that an increased LD can (1) prolong APD via the enhanced inflow of ChR2 current, and extend durations of twitch and internal sarcomere deformation through K⁺ currents and CaTnC complexes during ECC; and (2) extend durations of myocyte twitch and deformation and prolong APD through the inward ChR2 current and its inward rectification kinetics during MEF. An enhanced LI do not change APD and the duration of twitch in both ECC and MEF due to the ChR2 current's rectification feature. The optogenetic technology should not be viewed only as a tool to characterize cardiac electrophysiological properties. It also shows potential in the regulation of ECC and MEF in myocardium. Our study indicates that the inward ChR2 current and its inward rectification feature affect significantly the durations of AP and twitch in both ECC and MEF. The roles of optogenetic actuation on both ECC and MEF should be considered in future cardiac computational optogenetics at the tissue and organ scale.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under grant number 82060332, the Key Research and Development Program of Zhejiang Province under grant number 2020C03016, and the Hainan Provincial Natural Science Foundation of China under grant number 820RC625.

Conflict of interest

The authors declare no conflict of interest.