An anisotropic network model of inter-cellular coupling 1 for robust simulation of cardiac electrical excitation in 2 structurally complex and heterogeneous large tissue 3 models

7 Inter-cellular electrical coupling in the heart is a major determinant of excitation patterns in health 8 and disease. Cardiac arrhythmias, in particular rapid arrhythmias such as tachycardia and fibrillation, 9 are characterised by spatially complex and heterogeneous conduction patterns. In disease 10 conditions, these complex excitation patterns may result from electrophysiological and structural 11 remodelling, for example, remodelling of inter-cellular gap junctions and/or the proliferation of 12 fibrosis. Evaluation of the precise mechanisms by which local tissue structure determines global 13 arrhythmic excitation patterns is a major challenge, yet may be critically important for the 14 development of effective treatment strategies. Computational modelling is a viable tool to address 15 this challenge, but the established approaches for organ-scale simulations are unsuitable to capture 16 the impact of local conduction heterogeneities. In this study, we present a novel network model of 17 inter-cellular coupling for anisotropic organ-scale simulation of electrical dynamics that enables 18 cellular connections to be heterogeneously interrupted. 19 The presented approach is both simple and powerful. Intercellular coupling strength is weighted 20 based on local myocyte orientation and can be easily modified by sampling from continuous 21 distributions and/or removing individual cellular connections entirely; both approaches can be 22 differentially applied to axial and transverse cellular connections for full control over the conduction 23 substrate. 24 Preliminary simulations demonstrate the value of such an approach and indicate that perturbing 25 cellular connections in this manner can facilitate lower global conduction velocities and capture 26 wave breakdown and the development of re-entry in a way that is not possible with the previously 27 established approaches. Therefore, we have presented a useful model to simulate the impact of 28 local conduction heterogeneity in organ-scale models that expands the scope and possibilities for 29 computational models to elucidate arrhythmia mechanisms and generate genuine subject-30

The presented approach is both simple and powerful.Intercellular coupling strength is weighted based on local myocyte orientation and can be easily modified by sampling from continuous distributions and/or removing individual cellular connections entirely; both approaches can be differentially applied to axial and transverse cellular connections for full control over the conduction substrate.
Preliminary simulations demonstrate the value of such an approach and indicate that perturbing cellular connections in this manner can facilitate lower global conduction velocities and capture wave breakdown and the development of re-entry in a way that is not possible with the previously established approaches.Therefore, we have presented a useful model to simulate the impact of local conduction heterogeneity in organ-scale models that expands the scope and possibilities for computational models to elucidate arrhythmia mechanisms and generate genuine subjectspecificity, which may be particularly relevant in disease conditions.

Introduction
Cardiac arrhythmias are a major cause of the morbidity and mortality associated with cardiovascular disease (CVD).Loss of the regular rhythm of the heart, whether rapid and irregular (i.e., tachycardia or fibrillation) or abnormally slow (i.e., bradycardia), can substantially reduce cardiac output, drastically affecting the delivery of blood carrying oxygen and nutrients to the vital organs body, including the brain; acute arrhythmia events can be fatal and, indeed, sudden cardiac death accounts for a substantial proportion of mortalities associated with CVD (Virani et al., 2021).
Electrical and structural remodelling both occur over the progression of multiple CVDs, including heart failure and atrial fibrillation, and it is well established that this generally promotes both the initiation and sustenance of arrhythmia (Dhein et al., 2014;Han et al., 2021).Recent studies have highlighted the importance of structural remodelling, which may include hypertrophy, myocyte orientation disorganisation, reduced and heterogeneous cellular connections, and remodelling of the extracellular matrix including the proliferation of fibrosis, to both the triggers and substrate necessary for arrhythmia, contributing to conduction block and enabling sustained re-entrant circuits (e.g.Benoist et al., 2014;Prabhu et al., 2018).Experimental analysis of the local and global effects of structural remodelling is challenging, as the spatial scales required to image how local excitation depends on inter-cellular coupling and fibrosis preclude the global picture of whole-heart activation and arrhythmia patterns, in addition to the challenges of obtaining overlaid structure and function of the same tissue sample.
Multi-scale computational modelling has proved an invaluable tool in assessing the relative contributions of electrical and structural remodelling to arrhythmia, and enabling elucidation of the underlying mechanisms; such insight can help to drive novel pharmacological and surgical interventions to prevent, manage, or treat CVD in a diverse population (Lamata et al., 2014;Vadakkumpadan et al., 2009;Whittaker et al., 2019).Despite the successes of tissue-scale modelling of cardiac electrical dynamics (Bakir et al., 2022;Colman et al., 2022), the approaches used for organ-scale simulations are not well suited to capturing the impact of local heterogeneous conduction and conduction barriers associated with remodelling and, in particular, interstitial fibrosis.Given that these features promote arrhythmogenesis and substantially determine arrhythmia dynamics, this limits the capability of the models to understand the fundamental underlying mechanisms as well as generate genuine subject-specific models at this scale.
The two most commonly used approaches for simulation of organ-scale cardiac electrophysiology are the finite difference method (FDM) and finite element method (FEM), referring to different mathematical approaches to discretise the mono-or bi-domain reaction-diffusion equations which describe spatial diffusion of electrical or chemical activity (Benson et al., 2021).FDM is the simplest approach but has limits in regards to stability in highly heterogeneous or anisotropic conditions.FEM is more complex and robust but is also not well suited to modelling interrupted cellular connections.Moreover, neither model well reproduces rapid regular pacing induced arrhythmia (often, S1S2 pacing timed within a vulnerability window or cross-field stimuli are required to induced re-entry using these approaches) nor small, very slow conducting re-entrant circuits (Campos et al., 2019).These established methods are unsuitable for simulation of these features primarily because they are derived from mathematical discretisation of a spatially continuous system of equations, yet cardiac tissue is inherently discrete at a sufficiently large spatial scale (that of the myocyte; 10-100 µm) that this discreteness may be important for dynamics, in particular in disease conditions.

3
Recently, cellular-scale tissue models have been developed which describe much more accurately inter-cellular coupling (e.g., Jaeger et al., 2021).However, large-scale tissue models, e.g., of the whole atria or ventricles, are computationally intractable at the cellular scale: spatial resolutions are typically 0.2mm or larger in each direction (corresponding to of order 1000 myocytes), providing a feasible number of nodes or elements on which to solve the system.Such models are constructed from high resolution imaging of the heart, including, where possible, reconstruction of myocyte orientation anisotropy (which may be from diffusion tensor MRI, or image-analysis of e.g., CT scans).It is not clear how these cellular scale coupling models could be applied in the context of these reconstructions.There is therefore a strong motivation to develop an approach to describe intercellular coupling at this larger than cellular resolution scale which can be easily derived from anatomical/microstructure data that are readily available and used in current models, while also exploiting the inherent discreteness of cardiac tissue and enabling simulation of highly heterogeneous media, including local barriers to conduction.Such a model would offer the possibility to reveal new insight into the mechanisms of cardiac (dys)function and develop ever more accurate models of the heart, including for genuine subject specificity.
In this study, we develop an anisotropic network model to achieve this goal.This manuscript will discuss the derivation of this approach and perform some demonstration simulations to illustrate how it may be used to model heterogeneous conduction conditions, in particular, the reduction to cellular coupling associated with disease and in the presence of increased fibrosis.

Overall approach
Cardiac myocytes are electrically coupled through intercalated discs which are connections between cells that contain gap junctions enabling inter-cellular flow of ionic currents (Kléber and Jin, 2021).Myocytes are elongated in structure (with a length approximately five-to-ten times their diameter) and gap junctions are preferentially located at the longitudinal (or axial) connections, forming fibre-like strands of coupled myocytes.Coupling is also observed in the direction perpendicular to myocyte orientation but is less than in the axial direction, and the conduction velocity in the axial direction is typically at least twice that in the transverse direction (Dhein et al., 2014).This feature of differentiated axial and transverse inter-cellular coupling -and the ability to interrupt them individually -underlies the philosophy our proposed novel approach.
The approach operates on a structured gird of nodes i.e., a square (2D) or cuboid (3D) lattice.This represents the simplest discrete reconstruction of cardiac images, as opposed to an unstructured gird comprised of polygonal elements and vertices, and corresponds to the same system on which the FDM approach is applied.In this system, we present the inherently discrete equation governing the voltage differential for some node, i Where Igap (i,n) is the current associated with gap junction n that is connected to node i, and the sum over n corresponds to all junctions associated with node i (in 2D this is a sum of up to eight junctions; in 3D it is a sum of up to 26 junctions).It is common to absorb the membrane capacitance, Cm, into the formulation of the ionic currents, through the definition of the maximal conductances in units of nS/pF.The same will be applied here for the conductance of the gap junctions, and so equation (1) simplifies to: (2) The junctional current for some junction, n, between two adjacent nodes, i and j, is defined by: (3) Where ggap n is the conductance of gap junction n.The gap junctional current for the nodes i and j, the term that contributes to the sum in equations ( 1) and (2), is therefore given by: (4) (5) Note that we use the terms "gap junction", "junctional current" and "junctional conductance" broadly here, as they may not correspond to individual gap junctions between individual myocytes, but rather the connections between nodes.Moreover, this approach is not a formulation of gap junctional dynamics itself, and the conductance of the junction will be assumed to be constant (i.e., not voltage-or time-dependent).Rather, we present a method to discretise the tissue model in order to determine the magnitude of this gap junctional conductance between neighbouring nodes based on local myocyte orientation.
We must now determine how nodes are connected to form junctions, and how to derive the gap junctional conductance, ggap n based on local myocyte orientation.Throughout the remainder of the methods section for brevity, only the minimum number of equations necessary to understand the construction of the model will be shown, exploiting symmetries (e.g., x-directions will be shown but not the symmetric y and z equivalents); full equations are provided in the Online Supplement.
2.2 Derivation of the approach in 2D

Axial and transverse connections
In a 2D structured grid, each node can be connected to eight other nodes, corresponding to four along the axes (±x, ±y) and four along the diagonals (±xy ++ , ±xy +-) (Fig 1A).Due to the different directions along each of these lines being indifferent (+x is the same as -x), this reduces to four independent connection directions: x, y, xy ++ and xy +-.Note that xy ++ refers to the diagonal where the sign on x and y is identical (+x and +y, or -x and -y) and xy +-to the diagonal where the signs are different.
The main concept of this approach is that, for each node, the contribution to the connection to its neighbours in each of these directions will be determined from weighting terms defined by the local myocyte orientation.Defining ga as the conductance of a gap junction in the axial direction (along the primary myocyte orientation) and gt as the conductance of a gap junction in the transverse direction (i.e., orthogonal to the primary myocyte orientation), the nodal junctional conductance in each direction (gei node ) can be described as: And similarly for the y and xy +-directions (Online Supplement), where gxx node is the contribution of the node to the gap junctional conductances in the x-direction, and equivalently for all other directions, Wxx is the weight towards the x-axis for the axial component, Wxx t is the weight towards the x-axis for the transverse component (and equivalent for all other directions), and Dx and Dy refer to the discretisation space step in each dimension and are included from geometrical arguments.Due to considerations regarding the relationship between coupling strength and the spatial step (see Discussion), it may also be desirable to define these conductances independent of the spatial steps.In this case, it is important to retain the factor of 1/Ö2 in the diagonal terms, or otherwise account for the increased distance in the diagonal directions, for geometrical consistency.

Defining the junction conductance and currents
Before describing how the weighting terms are derived, it is important to relate the nodal conductances (equations 6-7) to the junctional conductance in equation (3).A junction is formed between all pairs of nodes which correspond to tissue, adjacent to each-other in any of the coupling directions.For some junction, n, the conductance, ggap n in equation ( 3), is given by the mean of the nodal conductances in the direction in which they are coupled (Fig 1B).For example, for two nodes, i and j, adjacent to each-other in the x-direction, forming a junction n ( 8) And equivalently for all other directions in which nodes may be coupled.There is no distinction here whether the connection is axial or transverse, so mixed connections (e.g., at an abrupt change in fibre) are possible.Every time a junction is created, maps must be created (for implementation purposes) defining which node is the positive and which is the negative contributor to this junction The gap junctional current for junction n, equation ( 3), is therefore: And these junctional currents are summed for each node as described in equations (2, 4 and 5).

Deriving the weights based on myocyte orientation
The weights that scale the contribution of the axial and transverse gap junctional conductances to each nodal directional conductance term (equations 6-7) must be defined.The myocyte orientation in 2D will always point in a segment between one axis and one diagonal direction (Fig 1C).Due to periodicities and symmetries in the trigonometric functions, we can consider and calculate the weight towards the axis and diagonal independent of in which segment the orientation is pointing.We can first calculate the angle from the x-axis, θ: where Ox is the x-component of the normalised orientation vector, and Oy is the y-component.This will always return an angle between 0 rad (when the orientation is along the x axis) and π/2 rad (when the orientation is along the y axis).We can then define a weight which linearly depends on this angle as a measure of where the orientation lies between the x or y axis (θ = 0 or π/2 rad, respectively) and the diagonal (θ = π/4 rad): i.e., Waxis is equal to 1 if the orientation points exactly along either the x or y axis (θ = 0 or π/2 rad, respectively), equal to 0 if it points exactly along the diagonal (θ = π/4 rad), and linearly varies between 0 and 1 based on the angle in-between (Fig 1C).The weight towards the diagonal is then simply: (14) A conditional algorithm could then be applied to determine to which axis and diagonal the weighting terms should be applied.Alternatively, and perhaps more concisely, these choices can be absorbed into a single general equation through the introduction of binary pointing parameters.Thus, we may define: where Px is a binary parameter which is equal to 1 if the orientation vector points between the x-axis and either diagonal and equals zero otherwise; Pxy+ is equal to 1 if the orientation points between the diagonal xy ++ and either axis.Note that Py and Pxycould also be introduced, defined as 1-Px and 1-Pxy+, respectively.From this, the weighting towards the four directions is then given by: This approach can be illustrated by considering some simple example cases (Fig 1D).

Transverse direction
A major feature of this approach is to apply the transverse coupling always at an orthogonal direction to the myocyte orientation (Fig 1C), as opposed to FDM, for example, where the basic isotropic coupling is always applied along the axes, independent of the myocyte orientation.In 2D the transverse weights are trivial to calculate: Which are the same as equations (17-19) with the pointing parameters swapped.

Derivation of the approach in 3D
2.3.1 Primary myocyte orientation weights in 3D In 3D each node can be connected to 26 neighbours (i.e., 13 independent directions given the +/directional indifference), corresponding to the axes and diagonals in each plane and the four corners (Fig 2A).The primary myocyte orientation vector therefore points in a segment which corresponds to a quadrant on the surface of a cube, i.e., contributing towards the weight of up to four directions at once (Fig 2B).Similar to the approach in 2D, we can consider how the weights towards these four directions depend on the angle(s), independent of which axes and diagonals define the quadrant.
The four directions (Fig 2B) correspond to the axis, the in-plane diagonal (dependent on θ1, defined in Fig 2B ), the elevation diagonal (dependent on θ2, again defined in Fig 2B ) and the corner (dependent on both θ1 and θ2).Note that both θ1 and θ2 are defined as the projection in the axis planes i.e., θ2 does not correspond to the spherical polar coordinate's elevation angle.We can then define the weight towards each of these directions as (Fig 2B ): where Waxis plane and Waxis elevation are defined from the two angles, θ1 and θ2, in the same way as the 2D model (dependent on which axes contribute to these calculations).The three angles (θxy, θxz, θyz) and corresponding weights (Waxis xy , Waxis xz , Waxis yz ) in each plane are defined by: and the pointing parameters in 3D become: (27) ) . 1−W axis elevation ( ) And similarly for the y-and z-directions (Online Supplement).Note that only one of Px, Py and Pz can be non-zero at any time (the orientation can only be pointing primarily towards one of these axes).
Pxy+, Pxz+ and Pyz+ can all be non-zero.In full, therefore, the weights towards the four quadrant directions are given by: Note that because only one of Px, Py and Pz can be non-zero at a time, these equations reduce immediately to the equivalent of equations (23-26).
The weighting towards each direction is defined by: And similarly for all other directions (Online Supplement).Note that the form of equation ( 36) results from the convention that when pointing primarily towards z (Pz =1), the xz diagonal is considered "in plane" and the yz diagonal considered the "elevation" direction, whereas when pointing towards x or y, the xz and yz are both considered the elevation direction.This is an arbitrary choice, and the results would be the same whichever choice is made.As with the 2D case, when the pointing parameters are factored into these equations, only a maximum of four terms will result in being nonzero: one axis, two diagonals, and one corner.

Transverse connections in 3D
In 3D, two transverse directions are necessary.These may correspond to the distinguishable "sheet" and "sheet normal" directions (LeGrice et al., 1995) if these data describing the laminar/sheet structure of cardiac tissue are available [e.g. from diffusion tensor MRI (Helm et al., 2005;Holmes et al., 2000;Hsu et al., 1998;Scollan et al., 1998;Tseng et al., 2003)] or to two indistinguishable (equal magnitude) transverse directions if only myocyte orientation data are available (i.e. in the absence of data describing the laminar/sheet architecture of cardiac tissue).The junctional conductance is therefore given by: (37) (39) And similarly for all other directions (Online Supplement).Similarly to the 2D case, if these are defined independent of the spatial step, the factors of 1/Ö2 and 1/Ö3 must be retained for the diagonal and corner directions, respectively.If differentiating between in-sheet and transverse to sheet, then gt1 > gt2, whereas if not imposing this distinction, gt1 = gt2 =gt and the equation(s) reduce to: (40 If implementation is based on three eigenvectors (e.g, from DTI), then the orientation components of the two transverse vectors directly determine Wei t1 and Wei t2 for each direction through the same equations (34-37) as the primary orientation vector (Fig 3).However, if only information on the primary eigenvector is defined, then we must define the transverse vectors from this (Online Supplement).

Implementing heterogeneous media
One main advantage, and indeed motivation for the development, of the above approach is to be able to modify cellular connections directly.For example, connexins (gap junctions) are heterogeneously expressed in tissue and may not exist between every adjacent myocyte.In disease in particular, highly heterogeneous media may be observed, corresponding to heterogeneity in connexin expression, or increased fibrosis.The latter may be associated with a loss of cellular connectivity and changes to the anisotropy ratio of conduction velocity, implying differential loss of transverse and axial connections.
Implementing these heterogeneities in this network model may be approached in different ways, depending on the application.For example, a map can be passed in that simply scales local (nodal) ga and gt, in the same way as one might do when using FDM.However, we can also scale the junctions directly, or even remove them.In this study we demonstrate these different approaches and their potential impact on function, although it should be clarified that we are not constructing physiologically validated parameter combinations for different health and disease states but, rather, demonstrating how this model may be applied in these contexts.
This approach readily facilitates differential removal of axial and transverse cellular connections.One can define a threshold for each type of connection that determines the probability of a connection being removed.It is important to identify if a connection is axial (the contribution from both nodes is from the primary orientation vector), transverse (the contribution from both nodes is from a transverse vector) or mixed (one node contributes an axial component and the other a transverse).Once complete, random numbers can be generated and combined with the threshold to determine which connections are removed.Alternatively, as a more sophisticated approach, junctional conductance could be scaled by continuous factors, for example, randomly sampled from various user-defined distributions or based on analysis of variability in experiment.

Cellular dynamics and model parameters
Example simulations presented in this study were performed using cellular electrophysiology described by the hybrid minimal model presented in Colman (2019).This is a simple model parameterised for the magnitudes of ionic currents and voltage range observed in cardiac electrical excitation.It is useful and desirable to be able to define the connection parameter conductances (ga and gt) from the diffusion coefficient, D, to enable easy implementation to replace FDM, or similar, where D has been defined to match conduction velocity.A conversion factor can be obtained through analysis of a simple case (Online Supplement), yielding: (41) (42) As will be shown in results, this conversion gives the same conduction velocities as the FDM method in these simple conditions.Note that this linear relationship means gt can be defined from the anisotropy ratio (AR) and ga, in the same way as often for D1 and D2: In this setup, the units of ga and gt are the same as those for the maximal conductances of the ion currents (nS/pF), as the membrane capacitance (Cm) in equation ( 1) has been absorbed into these maximal conductance terms (equation 2).However, this model does not require this to be the case and enables coupling between cells of different Cm while still conserving current; one could define ga as an absolute conductance (nS) and maintain inclusion of the factor Cm -1 when implementing equation (1).
g t = g a / AR 3 Results

Connection maps in 2D
Firstly, cellular connections defined by the approach in different myocyte orientations are examined to ensure that they correspond to the orientation vectors.For three simple cases in 2D (myocyte orientation points exactly in x, exactly in xy ++ , and in between xy ++ and y; Fig 4A ), it is clear that the cellular connections follow the orientation.It is worth noting that when the orientation points either along the axis or diagonal, the junctional conductance in that direction is exactly equal to ga (the weight is equal to one) and the transverse coupling is applied also only in one direction (perpendicular); in these special cases, only two directions have a non-zero weighting (one axial and one transverse; Fig 4A).However, when the orientation points between these directions, the maximum junctional conductance will be less than this, as it is weighted between two directions; the transverse weight is similarly split between two directions and, overall, the conductance is weighted between all directions.
Cellular connections are also shown for the more complex condition of having a spatially-varying myocyte orientation, and in the case where a set proportion of axial and transverse cellular connections have been removed (Fig 4B).It is worth explaining the patterns that are observed in cellular connections, for clarity of interpretation of how the model works.If the orientation points exactly in x, gxx node will be exactly equal to ga (Wxx = 1; Wxx t = 0), gyy node will be exactly equal to gt (Wyy = 0; Wyy t = 1) and there will be no diagonal connections (Wxy++ = Wxy+-= 0).As the orientation is rotated towards the diagonal, the axial weight for the x direction will linearly reduce from 1 to 0 (while the weight for the diagonal linearly increases from 0 to 1) and the transverse weight for the y direction also linearly reduces from 1 to 0. Thus, exactly at the diagonal, gxx node and gyy node are 0 for both axial and transverse.As the orientation is rotated further towards y, the transverse weight for x increases from 0 to 1 along with the axial weight for y.Thus, gxx node is smallest when the orientation points towards the diagonal, increases to gt as the orientation rotates towards y, and increases to ga as the orientation rotates towards x.

Connection maps in 3D
Connection maps in each direction are shown for two 3D geometrical models: a human ventricular wedge (Benson et al., 2007) where only the myocyte orientation was given (and so the two transverse directions are calculated within the model; Online Supplement), and a full reconstruction of rat bi-ventricular geometry (Whittaker et al., 2019) wherein all three orientations describing cardiac tissue structure (myocyte, sheet and sheet normal, corresponding to the three DT-MRI eigenvectors) were given (Fig. 5).

Verification of conduction patterns and comparison to FDM
The network model was parameterised to match FDM implementations in regards to axial and transverse conduction velocity (see "Methods: Cellular dynamics and model parameters").Conduction patterns in 2D and 3D between the network model and standard FDM implementations are compared (Fig 6).There are some notable differences between the two implementations, but the overall patterns and activation time match well, and it is not necessarily the case that FDM represents a "ground truth".Note that the spatial pattern in both methods when the orientation is diagonal is not a symmetrically rotated version of the pattern when the orientation is in-axis.This is discussed further in the Online Supplement.

Conduction patterns in heterogeneous media
The impact of removing cellular connections at defined probability thresholds is demonstrated using a 2D model, implementing the heterogeneous vector field illustrated in Fig. 4, across total spatial extents of 50´50 nodes and 200´200 nodes, enabling direct relation to visualised connection maps (Fig. 7 and Supplementary Videos 1 and 2).It is clear that removing axial and transverse connections causes different spatial patterns to emerge, and that the activation patterns correspond to the differential loss of connections in these directions.Non-uniformity in the propagating wavefront is observed in all conditions where cellular connections are removed, and arguably more closely represents experimental data (e.g.optical mapping) than the idealised control situation.Patterns in media with sampled distributions of scale factors are shown in the Online Supplement.

Wave breakdown and re-entry in heterogeneous conduction media and low intercellular coupling regimes
Preliminary simulations implementing heterogeneously removed cellular connections reveal that these conditions are sufficient to promote propagation wave breakup, which, in combination with rapid pacing, can lead to the development of re-entry which is not observed with a parametermatched global reduction to inter-cellular coupling in FDM (Fig. 8 and Supplementary Videos 3 and 4).It should be clarified that the regimes shown in these cases correspond to substantial loss of connections and thus slow global conduction velocities or activation times, and this is reflected in square-like waves in the FDM implementation.

Summary
We have presented a novel approach to construct a network model to describe cardiac inter-cellular electrical coupling based on structured-grid and local myocyte orientation data.The approach works by calculating coupling weights between adjacent nodes based on the local myocyte orientation.We have demonstrated that this approach leads to expected conduction patterns in idealised and non-idealised tissue.We have further demonstrated that this approach readily facilitates the direct manipulation of cellular connections, enabling complex conduction patterns to be simulated, such as observed with high levels of fibrosis.The model conserves current, is highly stable to complex media, and naturally accounts for boundary conditions.

Value of the model
The model presented has some advantages compared to the commonly-used FDM and FEM in regards to computational efficiency, simplicity, numerical stability, ability to solve with simple numerical methods, and mass conservation.However, the main advantages -and primary motivation for the development of the model -are: (i) the ability to directly modify or remove cellular connections, resulting in representation of barriers to conduction between excitable regions that may more accurately describe fibrosis and heterogeneous connexin expression and permits robust globally slow conduction.
(ii) the setup readily enables direct integration with models describing gap-junctional conductance dynamics as a function of voltage, which recent studies have shown to be important for functional cellular coupling (Bragard et al., 2021;Hawks et al., 2019).
Example simulations demonstrate that removing cellular connections can produce complex, heterogeneous and highly anisotropic conduction patterns which may degenerate into arrhythmic dynamics, and that such a process can induce the same global effect of reduced conduction velocity but with a distinct underlying pattern to a simple global reduction in cellular coupling.Due to the importance of electrotonic interactions on the generation of behaviour such as alternans (Colman, 2020;Huang et al., 2020) and afterdepolarisations (Alexander et al., 2022;Colman, 2019), it may be critically important in disease models to accurately capture the underlying coupling substrate, rather than implement a simple global reduction to cellular coupling.This is reflected in conduction patterns in low coupling regimes, where the FDM method, for example, produces unphysiological square-like waves.
Due to these features, we argue that this model presents advantages over FDM implementations in all conditions, and over FEM implementations for the specific context of heterogeneous cellular coupling.Implementation of the model readily enables direct replacement of previous implementations of the FDM (as the same geometrical and myocyte orientation reconstructions can used, and conduction parameters easily matched).For this purpose, the implementation is made available open-source in two different packages: i) packaged with the Multi-scale cardiac simulation framework (Colman, 2019) and ii) in a simplified, single-file implementation in C++ to be extracted and used in any context.Tools to generate heterogeneous connection maps are also provided with the code, all available at http://physicsoftheheart.com/downloads.htmland https://github.com/michaelcolman/

Comparison to other non-FDM or FEM approaches
A number of recent studies have presented alternatives to the established approaches for discretising cardiac tissue and solving the propagation of electrical excitation.However, as far as we are aware, no previous study has attempted to develop an anisotropic network model in 3D based on cellular connections at larger then cellular-scale resolutions in the same way as the present study.Network models of electrical coupling have been previously presented (e.g.Kazbanov et al., 2016;Ten Tusscher and Panfilov, 2007;Vandersickel et al., 2018) but these models were not designed to implement complex anisotropy as presented in this model.Moreover, fibrosis and heterogeneous tissue were implemented in the same way as one might for FDM or FEM -that is, by removing excitable nodes in the tissue and/or globally reducing the diffusion coefficient.This is fundamentally different to our approach, wherein connections can be removed while maintaining each tissue node as excitable.We note that at a discretisation of 0.2 -0.35 mm, as is common in large tissue models, a single node corresponds to hundreds or thousands of myocytes and thus removing whole nodes from being excitable tissue may not be an accurate representation of fibrosis or barriers to conduction, whereas our preliminary simulations demonstrate that substantially different behaviour emerges in this approach compared to a global reduction in diffusion.Jaeger et al., (2021) recently presented a sophisticated model of inter-cellular coupling at the individual myocyte scale, which accounts for the three domains of the extra-cellular, intra-cellular and membrane spaces.Such an approach enabled simulation of "micro-re-entry", which cannot be reproduced using larger scale discretisations.Other approaches, also at the cellular scale, have explicitly described the dynamical conductance of gap junctions (Bragard et al., 2021;Hawks et al., 2019;Hurtado et al., 2020) and demonstrated dynamic, voltage-dependent properties to be an important regulator of electrical conduction.These models, however, do not directly compare with the novel model presented in this study regarding motivation, spatial resolution, or implementation; rather, the approaches are complementary.A major advantage of our approach is that gap junctional conductance, or coupling strength between pairs of nodes, is an explicit parameter of our model; it would be trivial to modify this term with appropriate gating variables to describe the voltage-dependence of these junctions, or by expanding it into the tridomain description.We have presented an algorithm which can be applied to larger-than-cellular-resolution tissue models based on local orientation information by weighting contributions between pairs of nodes, and thus our approach can bridge the gap between these cellular-scale approaches and organ-scale simulations.

Limitations and further work
Whereas this method is based on "cellular connections", the connection between nodes in a discretised geometry does not directly correspond to individual myocyte connections, as each node/voxel occupies a volume that contains multiple myocytes.This limitation is equally true of FDM and FEM and any tissue implementation that is not at the direct cellular scale.This is because it is not currently feasible to solve large tissue models (including whole-chamber or whole-heart) while accounting for every individual myocyte, due to both computational intractability and resolution limits on whole-heart imaging data.However, this does imply that modifying "cellular connections" at this scale is not directly translatable to heterogeneous connections at the individual myocyte level.Nevertheless, it still enables implementation of functional barriers to conduction without requiring excitable nodes to be removed which we argue is a more accurate representation of the underlying substrate.To further bridge this gap, future studies could develop high-throughput image analysis approaches to derive this model from local, high resolution histological data (e.g., describing connexin expression and fibrosis) further to the whole-heart data for which it has been designed, and appropriate approximations could be developed in simulations to match models at different scales.
A further limitation, albeit another not specific to this model but rather of most approaches to modelling spatial coupling, is that the absolute coupling strength between two nodes is dependent on the spatial resolution in order to maintain a set conduction velocity (the spatial step size appears in equations 38-40 and 49-50).However, electrotonic coupling between myocytes is clearly not dependent on the spatial resolution of a geometrical model, and this may have important implications as changing the electrotonic coupling will affect cellular action potential dynamics, such as upstroke and the emergence of after depolarisations.Whereas not implemented in this initial model development, a potential solution to this problem could be to have an absolute coupling strength (independent of the spatial discretisation step), and introduce a delay in the voltage "seen" by neighbours such that the conduction velocity can be maintained without actually affecting electrotonic coupling strength.This will be considered in future work.
Finally, there is a choice of how to treat mixed connections in the model (where one node contributes an axial but the other a transverse connection), and it is not clear which approach would be the most accurate.For example, a mixed connection would be observed at an abrupt change in fibre orientation, where the axial conductance of one node is coupled to the transverse of its neighbour.As presented in this study, the junctional conductance was just an average of the two contributing conductances, independent of which type of connection they were.One could just as easily make a choice to model such junctions as either fully axial (axial connection only contributes to the junctional conductance; perhaps this corresponds to an abrupt change representing myocyte branching at smaller scales?) or fully transverse (perhaps representing one "fibre" ending at a junction with a perpendicular fibre).This is a choice which could be carefully explored in future studies.

Conclusion
We have presented a simple alternative formulation to the reaction-diffusion equation, which sums over gap junctional currents, rather than discretising the continuous gradients.Whereas the model superficially has a lot of similarities to simple FDM implementations, it also contains a fundamentally different underlying approach which presents a number of advantages, including model stability and conservation of current, but most notably in the ability to directly control cellular connections, and suitability for integration with dynamic, voltage-dependent gap-junction models.Both of these features enable more accurate and robust simulations of cardiac tissue to be performed, particularly in the disease state where tissue can become highly heterogeneous and anisotropic, opening new avenues for the elucidation of underlying arrhythmia mechanisms in disease.

Acknowledgements
The authors would like to thank Dr Dominic Whittaker for helpful conversations during the early development of the approach, and Dr Enric Alvarez Lacalle and Prof. Arun Holden for providing useful feedback on early versions of the manuscript.

Figure 1 .
Figure 1.Illustration of the model connections in 2D.A -illustration of the different coupling directions and terminology used in this paper.B -illustration of how the nodal and junctional conductance parameters are related and contribute to the total gap junctional current, for three nodes i, j, and k and two junctions n and m.C -illustration of the relationship between myocyte orientation angle and the weighting terms towards the axis and diagonal; in this particular example, the myocyte orientation is between the positive x axis and the xy ++ diagonal.D -illustration of the different values the binary pointing parameters take dependent on in which segment the myocyte orientation points.

Figure 2 .
Figure 2. Illustration of coupling directions and orientation weights in 3D.A -The 13 different coupling directions (nodal connections) in a 3D structured grid.B -illustration of how the weight towards each direction in a quadrant is calculated.

Figure 3 .
Figure3.Illustration of the transverse orientation vectors and the directions to which they contribute in 3D.Left -illustration of the three orientation vectors as given by imaging data; each vector points in a quadrant, contributing to up to four directions.Right -illustration of the definition of the two transverse vectors if only the primary eigenvector is given.From the axial orientation, a π/2 rotation is applied towards the z-axis to define transverse vector 1.The cross-product of these two vectors then defines transverse vector 2, which will always point in the x-y plane.Note that only contributions in the positive direction are shown for clarity: all coupling is in both directions along each vector, as illustrated in Fig 2B. g

Figure 4 .
Figure 4. Connection maps in 2D.A -Illustration of cellular connection conductance (indicated by the colour of the markers) associated with different (global) myocyte orientation directions, illustrated on a 10x10 grid for clarity in visualisation.B -Illustration of cellular connections associated with spatially heterogeneous anisotropy, demonstrating the myocyte orientation (left)and the cellular connections derived from this orientation, in both control (ga and gt are spatially homogeneous; upper) and in a case where 20% of axial and 80% of transverse connections have been randomly removed (lower).Note that, for clarity of visualisation, the junctional conductances defined in equations (38-40) have not been scaled by the 1/Dx etc factors, in order to normalise between axis and diagonal conductances.In the working model, for example, the junctional conductance in the diagonal direction is a factor of 1/Ö2 smaller than in the axis.

Figure 5 .
Figure 5. Illustration of directional conductances in a full bi-ventricular model, with all three eigenvectors given.Primary myocyte orientation streamlines are shown to provide context, with colour corresponding to either the x-or y-component (dependent on the view), along with the magnitude of the connection for each direction.Orientation and conductances are visualised on two slices, rather than throughout the volume, for clarity.The anisotropy ratio is 8:2:1 (i.e., ga = 4gt1 = 8gt2).As with Fig 4, for clarity of visualisation, the junctional conductances defined in equations (38-40) are not scaled by the 1/Dx etc factors, in order to normalise between axis and diagonal conductances.

Figure 6 .
Figure 6.Comparison between activation patterns in parameter matched FDM and network model simulations in 2D with different global orientations (fully in x, exactly in xy ++ ), and in the heterogeneous vector field shown in Fig 4, as well as in the 3D ventricular wedge model.Simulations with global orientation direction were performed on a 300´300 grid at a spatial resolution of Dx = 0.25 mm with coupling parameters of D1 = 0.4 mm/ms and D2 = 0.1 mm/ms (FDM method) and ga = 1.6 nS/pF and gt = 0.4 nS/pF (network model).Simulations using the vector field for myocyte orientation were performed on a 200´200 grid at a spatial resolution of Dx = 0.2 mm with coupling parameters of D1 = 0.2 mm/ms and D2 = 0.05 mm/ms (FDM method) and ga = 0.8 nS/pF and gt = 0.2 nS/pF (network model).Simulations on the 3D human ventricular wedge were performed at a spatial resolution of Dx = Dy = 0.2125 mm and Dz = 0.25 mm with coupling parameters of D1 = 0.1171 mm/ms and D2 = 0.0130 mm/ms (FDM method) and ga = 0.55 nS/pF and gt = 0.061 nS/pF (network model).

Figure 7 .
Figure7.Conduction patterns in different connection-removal conditions.A -Connection maps and activation patterns in a 50´50 grid, where the pattern can be somewhat directly visually related to the connection maps.Control refers to the condition where all connections are present, and the bottom row is the condition where 20% of the axial and 80% of the transverse connections are randomly removed.B -connection maps and conduction patterns in four different conditions (control and cases where 50%-50%, 20%-80% and 80%-20% axial-transverse connections have been removed) in a 200´200 grid, where direct visual relation to the connection map is more difficult, but the tissue has sufficient area to permit a reasonable activation time.

Figure 8 .
Figure 8. Illustration of wave-break and re-entry occurring in the network model with reduced connections.A -Voltage snapshots and activation pattern in the 2D model with the same vector orientation field as previous figures, paced from one edge of the tissue.B -Voltage snapshots in two different pacing scenarios (edge pacing for 5 stimuli, upper, and centre focal pacing for 2 stimuli, lower) on the same model as A. Sites of conduction block and wave propagation direction are labelled for clarity.
This work was supported by a Medical Research Council, UK, Strategic Skills Fellowship (grant no.MR/M014967/1) and Career Development Award (grant no.MR/V010050/1) awarded to MC, and a British Heart Foundation Project Grant (grant no.PG/16/74/32374) awarded to APB.