Theoretical Aspects of Resting-State Cardiomyocyte Communication for Multi-Nodal Nano-Actuator Pacemakers

The heart consists of billions of cardiac muscle cells—cardiomyocytes—that work in a coordinated fashion to supply oxygen and nutrients to the body. Inter-connected specialized cardiomyocytes form signaling channels through which the electrical signals are propagated throughout the heart, controlling the heart’s beat to beat function of the other cardiac cells. In this paper, we study to what extent it is possible to use ordinary cardiomyocytes as communication channels between components of a recently proposed multi-nodal leadless pacemaker, to transmit data encoded in subthreshold membrane potentials. We analyze signal propagation in the cardiac infrastructure considering noise in the communication channel by performing numerical simulations based on the Luo-Rudy computational model. The Luo-Rudy model is an action potential model but describes the potential changes with time including membrane potential and action potential stages, separated by the thresholding mechanism. Demonstrating system performance, we show that cardiomyocytes can be used to establish an artificial communication system where data are reliably transmitted between 10 s of cells. The proposed subthreshold cardiac communication lays the foundation for a new intra-cardiac communication technique.


Introduction
The heart's function is dependent on cardiomyocytes contracting in a coordinated fashion when electrically stimulated by the conduction system ( Figure 1a). The electrical activity starts at the Sinoatrial (SA) node-a node of specialized cardiomyocytes that initiates a synchronized electrical impulse. The SA node is a natural pacemaker and the electrical activity spreads to the right and left atria, depolarizing them to contract. The impulse spreads to the ventricles through the Atrioventricular (AV) node, the right bundle branches (RBBs) and left bundle branches (LBBs), and the Purkinje fibers. Electrocardiogram (ECG) is used to record cardiac electrical activity as a combination of all action potentials generated by the nodes and the cardiomyocytes.
In the presence of heart muscle damage, heart conduction may be disturbed, and artificial pacemakers are needed to re-establish regular cardiac operation [1,2]. In our recent paper, we discussed state-of-the-art pacemakers, and proposed a conceptual nano-actuator-network-based leadless pacemaker to overcome limitations imposed by battery longevity [1] (Figure 1b). Such a leadless device would pace numerous parts of the heart, using nano-actuators inter-connecting with individual cardiomyocytes, perform basic stimulation tasks by injecting current to the cytosol, and work in synchrony to optimize the energy used by individual batteries in the devices. Evoked electrical impulses/action potentials from actuated cardiomyocytes could then coordinate contraction throughout the heart muscle and lead to a normal heartbeat.  In this work, we explore the latter approach and use of membrane potential perturbations generated and propagated by cardiomyocytes when external stimuli and/or ion exchange occurs between the intra-and extra-cellular environments. In the considered scenario, if a stimulus applied to a targeted cardiomyocyte by a nano-actuator paces the heart, the cardiomyocyte should respond with action potentials. The action potential occurs when the membrane potential reaches the specified threshold values (typically −60 mV). Conversely, if a stimulus applied to the targeted cardiomyocyte by the nano-actuator is used to transmit data to another nano-actuator(s), the cardiomyocyte should respond by creating subthreshold membrane potentials. The subthreshold membrane potential thus refers to the membrane potential whose peak amplitude is below the specified action potential threshold. We envision the utilization of subthreshold membrane potentials in the time period between consecutive action potentials and use them as encoding signals ( Figure 2). In other words, transmission happens during the ventricular diastole phase. Of note, the theoretical framework presented in the following applies to different types of cardiomyocytes, e.g., those originating from the ventricles or atria. Nonetheless, by adopting a set of cell-specific parameters for the numerical simulations, we present numerical results that applies to the ventricular cardiomyocytes only.  Cardiac infrastructure indeed developed during evolution to conduct the signaling messages among cardiomyocytes and coordinate heartbeats. By utilizing encoding of subthreshold membrane potentials, the cardiac signaling system may be transformed into a more advanced cardiac communication system. We refer to the proposed communications paradigm between nano-actuators as Subthreshold Cardiac Communication.
In the work presented, we assume that the encoding subthreshold membrane potentials does not interfere with action potentials, nor affect normal heart function.
The rest of the paper is organized as follows. In Section 2, we propose a basic cardiac communication system and to linearize the cardiac cell membrane circuit by using the quasi-active method. In Section 3, we characterize the impacts of various noise sources on system performance. In Section 4, we provide numerical results. Finally, in Section 5, we discuss and conclude the study.

Subthreshold Cardiac Communication System
We considered the communication system between a transmitting nano-actuator and receiving nano-actuators within the multi-nodal pacemaker network. A small membrane patch connected to the transmitting nano-actuator is located in a selected compartment of the emitting cardiomyocyte and is used for stimulation. The emitting cardiomyocyte can respond to the provided simulation patterns with action potentials or subthreshold membrane potentials, which are both distributed forward through the unidirectional propagation channel [20,21]. The signaling/communication channel consists of cardiomyocytes connected by specialized gap-junction-like channels [22]. Gap junctions can be observed as aggregations of single dynamic and multi-functional channels, called connexins, which play a complicated and essential role in the entire conduction system of the heart [23]. A small membrane patch connected to the receiving nano-actuator is located in a selected compartment of the receiving cardiomyocyte. The receiving cardiomyocyte responds either with action potentials or subthreshold membrane potentials to the propagated signals. In this scenario, we inspect the performance of the described system when the transmitting nano-actuator sends encoded data during the subthreshold time-period by generating a stimulus to the emitting cardiomyocyte. Table 1 defines symbols used throughout the paper. Table 1. Defined symbols used throughout the paper.

Parameter Description Unit
Resistivity of membrane kΩ· cm 2 z l Equivalent longitudinal impedance (per unit length) kΩ/cm Z l Equivalent longitudinal resistivity kΩ· cm Output noisy voltage PSD mV 2 /Hz We opt to use one-dimensional cable theory to analyze the propagation of unidirectional subthreshold membrane potentials along cylindrically shaped cardiomyocytes [24,25]. The one-dimensional cable theory is widely used in the literature to model excitable tissues, e.g., nerve axons and skeletal muscle fibers [25]. Although cardiomyocytes generally form a strand consisting of individual cells with irregular shapes, the theory could potentially lead to inaccurate numerical results. Nonetheless, as this work lays the foundation for a new concept of biological communication paradigm, we believe that applying the one-dimensional cable theory reduces the complexity in this initial study. Besides, we use the quasi-active method to linearize the membrane's active channel kinetics into phenomenological impedances when subthreshold membrane potentials have small variations around the holding potential (The value of the holding potential refers to a specific value used as the baseline to determine fluctuations of the membrane potential.) [26][27][28][29]. The phenomenological impedance can have positive or negative components, depending on the difference between the holding potential and the reversal potential of different ionic channels [26]. The linearization is exclusively valid in the considered subthreshold regime where non-linearities, encountered in the creation of action potentials, are not expected to occur. The linearization, however, ensures us to use tools for studying linear systems and analyze the behaviour of this complex biological system.

Transfer Impedance
According to the one-dimensional cable theory, the propagation channel is equivalent to the cable that consists of the membrane impedance (per unit length), the intracellular impedance, and the gap junction impedance.
We denote the membrane impedance as z m , as shown in Figure 3a, and define using the resistivity of the membrane Z m where a denotes the radius of the cardiomyocyte strand. Aiming to define Z m , we first inspect the transmembrane current components on a membrane patch (i.e., the currents that depend not only on membrane potentials but also on opening/closing of ionic channels, e.g., sodium, potassium, and calcium, and the currents that depend solely on membrane potentials), and then apply the quasi-active method [26,27,29]. The linearized circuit for a membrane patch in a specific holding potential in the subthreshold regime ([−84, −60] mV) is shown in Figure 3b. Detailed linearization method used to define the circuit in Figure 3b is shown in Appendix A.1. Following Figure 3b, where The intracellular impedance and the gap junction impedance are commonly referred to as the equivalent longitudinal impedance (per unit area) [30,31]. We denote the equivalent longitudinal impedance as z l (Figure 3b), and define using the equivalent longitudinal resistivity Z l [24] where the typical value of Z l ranges from 0.6 to 36.6 kΩ· cm [24]. The equivalent impedance of the overall propagation channel, hereinafter referred to as the transfer impedance, finally stems from the one-dimensional cable equation [29] Equation (5) characterizes the membrane potential dynamics in the frequency domain at different propagation distances, where γ(j f ) is the propagation constant expressed as [24,27] γ where S V is the surface-to-volume ratio of the cardiomyocyte. With the propagation distance set to infinity as the boundary condition, i.e., V(∞, j f ) = 0, we define the transfer impedance of the channel [29] Z

Noiseless Input-Output Relation
Without any loss of generality, we consider the unipolar non-return-to-zero (NRZ) line code as the stimulus applied to the emitting cardiomyocyte and data to be communicated from one nano-actuator to another. Aiming to ease the formulation of a complete analytical framework, we characterize the NRZ line code in the frequency domain defining its Power Spectral Density (PSD) as where A is the applied current amplitude denoting transmission of bit 1 of duration T s , f is the operating frequency, and δ is the Dirac delta function. A = 0 denotes transmission of bit 0. Referring to (7) and (8), the output voltage PSD in the receiving cardiomyocyte is defined as

Noise in the Subthreshold Cardiac Communication System
Field stimulation and direct stimulation are the two approaches used for the stimulation of cardiomyocytes. With field stimulation, the microelectrode is not directly fixed to the cell. The stimulation affects the membrane through the extracellular solution, which leads to the generation of membrane potential fluctuations [32]. With direct stimulation-an approach used by the recently proposed nano-actuators [1]-the microelectrode is attached to the cardiomyocyte directly. This approach, however, induces the environmental disturbance in the form of input-dependent noise [33,34], and the membrane-related noise [35]. We refer to the input-dependent noise as the encoding noise. Encoding-and membrane-related noise are denoted as N 1 and N 2 in the cardiomyocyte communication system as shown in Figure 4.

Message Source
Encoder +

Encoding Noise
With the effect of encoding noise reflected through i N 1 (t), the injected signal/current is where i Tx (t) denotes the noiseless component. i N 1 (t) has already been studied in the relevant literature and derived from an autoregressive random process w(t) as [34] i where denotes convolution.
In a complex z-domain, w(t) is discretized as W(z) = a 0 + a 1 z −1 + a 2 z −2 + a 3 z −3 + · · · + a n−1 z −n+1 , where a 0 , a 1 , a 2 . . . a (n−1) are the coefficients of the n-th order autoregressive random process. Considering W(z) as a linear filter and replacing z with e j2π f , we define the PSD of the encoding noise N 1 Finally, the PSD of the input affected by the encoding noise is

Membrane-Related Noise
Compared to encoding noise, the membrane-related noise is more complex and composed of the (1) voltage-gated channel noise induced by stochastic opening/closing of the voltage-gated channels, (2) shot noise induced by random ionic release, and (3) thermal noise induced by intrinsic circuit dynamics [36,37].
We opt to describe membrane-related noise following the rationale presented in [37][38][39], where the subthreshold neuronal membrane potential noise was characterized. Owing to the similar excitable properties of neurons and cardiomyocytes, we represent the membrane-related noise source as an equivalent Gaussian current source and denote with I N2 u in Figure 5. In the following, we characterize membrane-related noise per unit distance assuming that ionic channels are homogeneously spread over the cellular membrane.

Voltage-Gated Channel Noise
We model the effect of stochastic opening/closing of voltage-gated channels in unit distance via conductance variations by defining a noise current i 1 n 2 in µA/cm whereĩ K (t),ĩ Na (t) andĩ Ca (t) are the noisy currents produced by the random opening/closing of potassium, sodium and calcium channels, respectively.g K (t),g Na (t) andg Ca (t) denote conductance variations of respecting ionic channels around their steady-state values when the holding potential is V 0 m . E K , E Na and E si are the reversal potential of potassium, sodium and calcium, respectively. The conductance variation is a collective phenomenon of multiple channels, not a single channel, and depends on the length of the propagation channel and ionic channel density [40,41]. Due to the tiresome mathematical derivation, we derive the PSD of the corresponding ionic current components The PSD of the voltage-gated channel noise is theñ

Shot Noise
The shot noise is affected by the random ionic release when different types of cations depolarize the cellular membrane and generate the membrane fluctuation while propagating to the cytosol. We adapt the PSD of the shot noise from the Schottky's formula and define [36] where q K , q Na and q Ca are the charges of the moving potassium, sodium and calcium particles, respectively.Ĩ K ,Ĩ Na , andĨ Ca are the sodium, potassium, and calcium currents in the frequency domain, respectively.

Thermal Noise
Thermal noise is known as thermal agitation which stems from the random movement of electrical charges in the electrical systems. Thermal noise has a significant impact on the performance of the receiving cardiomyocyte [42]. The considered cardiomyocyte communication system suggested by us consists of different passive components, i.e., cytosol-related resistors generated by intracellular dynamics, and membrane-related resistors and capacitors generated by the phospholipid bilayer cell-membrane. However, as described in [37], we ignore thermal noise evoked by the cytosol-related resistors, and only consider membrane-related passive components to calculate the PSD of the thermal noiseS where k is Boltzmann constant and T is the absolute temperature. Ultimately, the overall effect of the membrane-related noise is now described as or, alternatively,S where Z(y, j f ) is given in (7) andS N2 u (j f ) in (18).

Noisy Input-Output Relation
For the linear cardiomyocyte communication system, the PSD of the signal received at the receiving cardiomyocyteS Rx (x, f ) is whereS Tx (j f ) is given in (13),S N 2 (x, j f ) in (19) and Z(x, j f ) in (7).

Numerical Simulations
Analysis of the performance of the subthreshold cardiac communication system relies on the computational action potential models. A general and uniform action potential model for cardiomyocytes generally does not exist because the ionic current components vary for different cell types. The verification of the model depends on the experimental data. In this study, we linearize the membrane into the primary circuit and study the subthreshold cardiomyocyte communication by using the Luo-Rudy model, which is based on the Hodgkin-Huxley-type formalism [43,44]. The Luo-Rudy model is one of the commonly used ventricular cardiomyocyte action potential models and considers the most typical ionic current components, particularly, the sodium channel function operating in the subthreshold regime. The Luo-Rudy model provides the coefficients of different ionic channels. We list the general parameters used in our simulation framework in Table 2. As of note, the channel density and conductance of single sodium and calcium ionic channels are taken from the Reuter et al.'s experimental work [45][46][47], whereas the channel density and conductance of single potassium ionic channels are taken from the Shibasaki's experimental work [48]. Other parameters mainly originate from the Luo-Rudy model related works [27,43].
The linearization depends on the holding potential, which can be any value in the subthreshold range. In the simulation framework, we set the holding potential to be equal to the resting potential of −84 mV since it enables a broader amplitude range for data transmission: the stimulation range is 24 mV when the threshold is −60 mV. Table 3 lists the circuit parameters obtained by applying the method in Appendix A.1 to linearize the membrane at the selected holding potential; the parameters used as an input to the linearization method are from Table 2. The linearized values of phenomenological resistors and inductors of sodium, potassium, and calcium are negative as the holding potential is smaller than their reversal potentials [26].  Figure 6 shows changes in the transfer impedance, that we compute from (7), against the system frequency and propagation distance between the transmitting nano-actuator and the receiving nano-actuator. As intuitively expected, the transfer impedance decreases with an increase in both frequency and distance. According to the peak transfer impedance value and given action potential threshold, we determine the maximum stimulation current of 3.81 nA which can be applied to the cells in the data transmission mode of the cardiac system. Selecting proper transmission rates in time bins intended for data transmission (Figure 2) is important to minimize intersymbol interference at the receiver point. For the system performance demonstration, we select the bit transmission rate of 5 bit/s to plot eye diagrams in Figure 7 considering different system parameters. An eye diagram is a tool used in communications engineering for the evaluation of the combined effects of inter-symbol interference and channel noise on the performance of a baseband signal-transmission system. An open eye diagram corresponds to minimal signal distortion. A closed eye diagram corresponds to signal distortion. We plot eye diagrams at different distances considering the noiseless and noisy scenarios, respectively. We observe that in short-distance transmission systems the eye openings are wide with plentiful margin decisions at the receiver regardless of the effect of noise (Figure 7a-d). As expected, however, the effect of noise combined with the increased distance among nano-actuators appears as the closure of the eye diagram with low-amplitude signals (Figure 7f,h). In such scenarios, highly sensitive receivers are required by the system to measure the low-amplitude signals. Of note, present high-performance microelectrodes could ideally measure potentials as small as 0.015 mV [49]. The completely closed eye appears in a noisy considered scenario where twelve cells compose the communication channel (Figure 7h).
We now demonstrate coded data transmission over the channel composed of ten cardiomyocytes using Amplitude-Shift Keying. This modulation technique implies that binary 1 is represented by transmitting a fixed-amplitude wave for a fixed time duration; otherwise, binary 0 is represented. We select n = 2 bits to represent a symbol, meaning that M = 2 n = 4 different symbols could be encoded, i.e., symbols 'a', 'b', 'c' and 'd' represented by 00, 01, 10 and 11, respectively. We randomly generate 1000 bits to ensure that the symbols are equally represented. The probabilities of the symbols 'a', 'b', 'c' and 'd' are 0.240, 0.254, 0.250 and 0.256, respectively. Figure 8a,b show a small portion of the transmission bit stream consisting of 10 bits. Figure 8c shows a small portion of the received bit stream consisting of 10 bits without and with the consideration of noise sources, respectively. We observe the positive effect of the noise when bit 1 is transmitted. This effect stems from the input-dependent noise that apparently enhances the cardiomyocyte stimulation. Conversely, we observe the negative effect of the noise when bit 0 is transmitted. This effect stems from the channel noise that decreases the margin decision. Nonetheless, with the proper threshold selected according to the eye diagram, the receiver can successfully decode the received signals to corresponding 1/0. As shown in Figure 8d, the receiver decodes the signal successfully in both scenarios. As of note, the performance is highly dependent on the stimulation amplitude. Based on randomly generated and transmitted 10,000 bits, we evaluate the bit error rate (BER) when changing the stimulation amplitude starting from 1.5 nA, as shown in Figure 9. As expected, the BER decreases with the stimulation amplitude and reaches the minimum value of 5 × 10 −3 when the stimulation amplitude is 3.5 nA.

Concluding Remarks
We considered the communication system between transmitting and receiving nano-actuators within a multi-nodal pacemaker network. The subthreshold cardiac communication paradigm considered in this paper offers a potentially groundbreaking method for data transmission within the heart. The demonstrated transmissions showed that data could be successfully transmitted in the subthreshold domain over tens of cells only. The results are still insightful and provide initial information on how to distribute and deploy the relay nano-actuating node(s) in the multi-nodal pacemaker network. Combining the subthreshold cardiac communication system with the optimal stimulation methods may provide an energy-efficient pacing of cardiomyocytes.
The time bins when transmission can happen correspond to the duration of the ventricular diastole phase which is approximately 430 ms [50]. Based on the presented results and analyzed bit rates, a very limited amount of data could be transmitted. Nonetheless, for the essential function of a multi-nodal leadless pacemaker, where the nano-actuators primarily sense membrane potentials of the corresponding cells and assist in pacing, the proposed communication system could enable transmission of a status of the node's stimulation activity ensuring coordinated operation within the network.
The numerical results presented in the paper are based on the computational Luo-Rudy model of cardiac action potentials whose parameters are used to linearize the cardiac circuit. Though the Luo-Rudy model is not perfect, e.g., it does not consider the stochasticity of single ionic channels, it still serves as the basis for most computational models and studies involving myocytes and provides the coefficients for different ionic channels. For more precise results, in-vitro experiments are needed to obtain precise parameters for the specific cells. The experiments will also reveal the dynamics of the gap junction resistances as they change according to the potential between the gap junctions. In this study, we considered the resistance of the gap as constant.
Action potentials may affect the performance of the subthreshold transmission due to variations (jitter) of their initial times in different physiological environments. The action potential duration also varies in different physiological environments, which directly affects the length of the temporal bins intended for data transmission. On the contrary, prolonged data transmission in time bins between consecutive action potentials may affect action potentials. Besides, when multiple nano-actuators transmit data at the same time, the interference to the EGM (electrogram) and ECG may affect the performance of the proposed system. The EGM is used to measure the local signal in the tissue level. To measure the interference of encoding signals to the EGM, we would need to (1) apply more advanced 3D topological tissue models, (2) analyze the coupling between cells, (3) identify possible multiple paths between the transmitting and receiving nano-actuators, and, ultimately, (4) consider the position and timing/synchronization issue of nano-actuators. This indicates the direction for future work in this field.
Furthermore, future work should include the complex structure of the cardiomyocytes, such as syncytium structure or network structure, together with the timing of signal transmission between nano-actuators and the gateway/hub. Ultimately, in-vitro and in-vivo experiments are urgently needed to generate more precise circuit models and obtain real data on subthreshold membrane potentials propagation which will be used to verify the numerical results presented in this paper.

. Membrane Linearization
We linearize the sodium-, potassium-and calcium currents in cardiomyocytes following the same methodology presented in [25]. The sodium current is expressed as [43] whereḡ patch Na is the maximum conductance of sodium channels, v is the membrane potential, E Na is the reversal potential, η patch Na is the sodium channel density on unit area and it is 10 in the simulation framework, and γ Na is the conductance of single sodium channel. m, h and j represent the sodium active channel parameter, the sodium inactive channel parameter, and the sodium slow inactive channel parameter, respectively. They are functions of time and membrane potential and indicate the subunits of the sodium channel.
The membrane potential can be linearized around small changes denoted with δ. Accordingly, the sodium current variation can be expressed as where m v , h v and j v are the parameters at steady-state. In (A3), δI Na1 is constant becauseḡ Na , m v , h v and j v are constant in the steady state. δI Na2 , δI Na3 and δI Na4 indicate the small current changes of the subunit m, h, j, respectively. We consider m variation in the following and extend results to h and j. The active subunit m has two states: open and closed. m changes with time as where α m and β m are the coefficients of the subunit from closing state to opening state and opening state to closing state, respectively. α m and β m are functions of v that change for a small variation as We then yield where p ≡ d/dt. Finally, we yield and b ≡ α m + β m . For δI Na1 , we yield δv δI Na1 which we abstract as a resistor generated by sodium channels For δI Na2 and according to (A8), we yield δv δI Na2 which we abstract as a resistor r m in serial connection with an inductor L m Similar to δI Na2 , we define resistor and inductor parameters for δI Na3 and δI Na4 .

Appendix A.2. Derivation of Current Noise PSD
Similar to neuronal ionic channels, the sodium, potassium, and calcium channels in cardiomyocytes can be considered as finite-state Markov chains. Considering the sodium channels, the autocovariance of the corresponding current is derived from [38,51,52] where η line Na = 2πaη patch Na is the sodium channel density in unit length, γ Na is the conductance of a single channel, and m v , h v , and j v are steady-state values when the membrane potential is v. P Na,0|0 (t) is the conditional probability of all the subunits of sodium channels to open at time t = 0 P Na,0|0 (t) = (m v + (1 − m v )e −t/τ m ) 3 where τ m , τ h and τ j are the time constant of m, h and j, respectively. By using the Wiener-Khinchine theorem, the current PSD is where A = η line Na γ 2 Na (v − E Na ) 2 m 3 v h v j v . Similar to sodium, the autocovariance of the calcium current is where η line Ca = 2πaη patch Ca is the calcium channel density in unit length, γ Ca is the conductance of a single channel, E si is the reversal potential of calcium, and d v and f v are steady-state values when the membrane potential is v. P Ca,0|0 (t) is the conditional probability of all the subunits of calcium channels to open at time t = 0 The autocovariance of the potassium current is where η line is the potassium channel density in unit length, γ K is the conductance of a single channel, E K is the reversal potential of potassium, and X v is steady-state values when the membrane potential is v. P K,0|0 (t) is the conditional probability of the potassium activation channels to open at time t = 0 P K,0|0 (t) = X v + (1 − X v )e −t/τ X , where τ X is the time constant of X. The current PSD is