Modeling electrodiffusion in cardiac intercalated disc nanostructures

Modeling electrodiffusion in cardiac intercalated disc nanostructures Ena Ivanovic1, Marina Vannuci1,2, Nicolae Moise3, Seth H. Weinberg3 and Jan P. Kucera1 1 Department of Physiology, University of Bern, Bern, Switzerland 2 Graduate School for Cellular and Biomedical Sciences, University of Bern, Bern, Switzerland 3 Department of Biomedical Engineering, The Ohio State University, Columbus, OH, USA This is the MATLAB source code permitting to replicate the results presented in our article entitled ”Electrodiffusion in cardiac intercalated disc nanostructures alters cell-cell action potential transmission via ephaptic coupling: a model study” that was accepted for publication in the Journal of Physiology“. Download the ZIP file “CodeAndMeshes.zip” and decompress all files into a dedicated folder. Make sure that all files are located in the same folder. Details regarding the operation of the scripts and functions are presented briefly below. More details are provided as comments at the beginning of the scripts/functions as well as throughout the code. The code was developed with the 2019a version of MATLAB and should also run with newer versions. Use in terms of the Creative Commons CC-BY license (Attribution 4.0 International). Main scripts with their auxiliary functions Main_MATLAB.m Main_COMSOL.m Main_OSU.m The results shown in Figures 1, 2, 4, 5 and 6 can be replicated using Main_MATLAB.m. To replicate Figure 3, use Main_COMSOL.m. To obtain the results with the tortuous mesh shown in Figure 7, use Main_OSU.m. Of note, the difference between these files resides mainly in the generation of the finite element mesh. Main_MATLAB.m generates the mesh within the code, while Main_COMSOL.m and Main_OSU.m load an external file (MeshCOMSOL.mat and MeshOSU.mat, respectively). The mesh data in MeshCOMSOL.mat were generated using COMSOL Multiphysics (version 6.3) while the mesh data in MeshOSU.mat were generated using dedicated code at Ohio State University (OSU). At the beginning of each script, a number of parameters and flags are defined (as MATLAB variables, see Table 1) and the user must adjust them accordingly as described below. Table 1: Variables (note that not all are available in the three different main scripts) Variable Possible Values Definition Default value NaC 0, 1 Na+ channel distribution: uniform (0) or clustered (1) 1 PATTERN 1, 2, 3, 4, 5, 6 Na+ channel cluster distribution pattern (in the flat intercalated disc) 1: Small central cluster 2: Large central cluster 3: 5 central clusters (3rd row of Figure 6) 4: 7 clusters (4rd row of Figure 6) 5: 13 clusters (5th row of Figure 6) 6: 3 groups of 5 clusters (6 th row of Figure 6) 1 CleftWidth_nm Between 10 and 100 Cleft width in nm outside the perinexus or gap junction plaque 30 gapfactor Between 0 and 1 fraction of gap junctional coupling relative to normal (e.g., 0.01 for 1%) 0 Dscaling 1 Scaling factor to simulate reduced diffusion due to extracellular proteins in the rim around Na+ channel clusters 1 DYNCONC 0, 1 Constant (0) or dynamically changing (1) ion concentrations 1 BLKOUTINA 0, 1 Outward Na+ current: active (0) or blocked (1) 0 GJp 0, 1 Gap junction distribution: uniform (0) or clustered (1) 0 PeriWidth_nm* Between 10 and CleftWidth_nm Perinexus width in nm 10 To reproduce the simulation results for each Figure, set these variables as follows: Figure 1: Open Main_MATLAB.m. For a uniform Na+ channel distribution (left column of Figure 1), set NaC=0. To regroup all Na+ channels in the centre of the intercalated disc (middle and right columns of Figure 1), set NaC=1. Choose the radius of the cluster by setting PATTERN=2 (middle column) or PATTERN=1 (right column, default). Leave the remaining variables at their default values (BLKOUTINA=0, Dscaling=1). To replicate the results shown in panel B, set CleftWidth_nm=40 and gapfactor=0. Obtain the black traces with DYNCONC=0 and the red traces with DYNCONC=1. At the end of the simulation, the delay between the upstrokes is shown in the MATLAB command window. To obtain the individual data points used to generate panel C, different simulations can be run with different values of CleftWidth_nm and gapfactor. Figure 2: Open Main_MATLAB.m. Set NaC=1 and PATTERN=1 to obtain a mesh with a small central Na+ channel cluster. Moreover, set Dscaling=1, GJp=0, gapfactor=0 and CleftWidth_nm=50. Similar to Figure 1, set DYNCONC=0 and DYNCONC=1 for results with constant and dynamic ion concentrations. To obtain simulation results where outward INa is blocked set BLKOUTINA=1. At the end of the simulation, the delay between the upstrokes is shown in the MATLAB command window. To obtain the individual data points used to generate panel C, different simulations can be run with different values of CleftWidth_nm and gapfactor. Figure 3: Open Main_COMSOL.m. To replicate the results in panel B set gapfactor=0 and CleftWidth_nm=50. Left column: DYNCONC=0 and BLKOUTINA=0. Middle column: DYNCONC=1 and BLKOUTINA=0. Right column: DYNCONC=1 and BLKOUTINA=1. Vary PeriWidth_nm between 10 and 50 nm in steps of 10 nm to reproduce the different curves. The simulations analyzed in panel C can be reproduced by setting Cleft_width_nm to values between 10 and 100 nm (in steps of 10 nm), PeriWidth_nm to values between 10 and Cleft_width_nm (in steps of 10 nm) and setting DYNCONC and BLKOUTINA accordingly. Figure 4: Open Main_MATLAB.m. Set NaC=1 and PATTERN = 1 to obtain a mesh with a small central Na+ channel cluster. To replicate the results in panel B set gapfactor=0, CleftWidth_nm=50 and vary Dscaling from 1 (default) to 0.2 in steps of 0.1. Set DYNCONC and BLKOUTINA as in Figure 3. The simulations analyzed in panel C can be reproduced by setting Cleft_width_nm to values between 10 and 100 nm (in steps of 10 nm), Dscaling to values between 0.2 and 1 (in steps of 0.1), and setting DYNCONC and BLKOUTINA accordingly. Figure 5: Open Main_MATLAB.m. Set NaC=1 and PATTERN=4 to obtain a mesh with multiple Na+ channel clusters. To replicate the results, set CleftWidth_nm=30, gapfactor=0.01, DYNCONC=1 BLKOUTINA=0 and vary Dscaling from 1 (default) to 0.2 in steps of 0.1. Figure 6: Open Main_MATLAB.m. Choose a Na+ channel distribution pattern (corresponding to the different rows of Figure 6) and set PATTERN accordingly (see Table 1). Set DYNCONC=1 and BLKOUTINA=0. To obtain the schematics shown in panel B, set gapfactor = 0 and to obtain the results in panel C, set gapfactor=0.01. Repeat the simulations for different values of CleftWidth_nm (from 10 to 100 nm in steps of 10 nm) and Dscaling (from 0.2 to 1 in steps of 0.1) to reconstruct the data. At the end of each simulation, the delay between the upstrokes is shown in the MATLAB command window. Figure 7: Open Main_MATLAB.m to replicate the results for a flat intercalated disc mesh with uniform Na+ channel distribution (NaC=0, Mesh 1). Open Main_OSU.m to replicate the results for the tortuous intercalated disc mesh, as follows. Mesh 2: set NaC=0, GJp=0 and Dscaling=1. Mesh 3: set NaC=1, GJp=0 and Dscaling=1. Mesh 4: set NaC=1, GJp=1 and Dscaling=1. Mesh 5: set NaC=1, GJp=1 and Dscaling=0.2, 0.3 or 0.4. To obtain the results depicted in panel B, set gapfactor=0.01, CleftWidth_nm = 20, DYNCONC=1 and BLKOUTINA=0. Then, execute the script. The program displays the mesh (and the boundary nodes) in a dedicated figure. The program also displays the membrane potential for the intercalated disc of cell 1 and cell 2, the sodium current densities for both cell membranes, the extracellular potential (in mV) in the cleft as well as ion concentrations in the cleft (ion species: sodium, potassium, calcium) as a colormap. At the end of the simulation, the program produces a figure similar to the ones presented in the article. Furthermore, the program outputs in the MATLAB command prompt the delay between the two upstrokes of the intracellular potentials. The script Main_OSU.m also generates movies, as detailed below. Movie S1: Open Main_OSU.m. Set the parameters as follows: NaC=0, GJp=0, Dscaling=1, gapfactor=0.01, CleftWidth_nm = 20, DYNCONC=1 and BLKOUTINA=0 (these are the parameters for Mesh 2). Execute the script, which by default will generate a movie file “Movie.mpg” in the working directory. Movie S2: Open Main_OSU.m. Set the parameters as follows: NaC=1, GJp=1 and Dscaling=0.2, gapfactor=0.01, CleftWidth_nm = 20, DYNCONC=1 and BLKOUTINA=0 (these are the parameters for Mesh 5). Execute the script, which by default will generate a movie file “Movie.mpg” in the working directory. The functions used by the three main scripts are briefly described below. For further details, see the comments in the corresponding MATLAB code. Functions related to the finite element method used in this work The partial differential equations were discretized in space using the finite element method with triangular elements and piecewise linear basis functions, resulting in systems of ordinary differential equations. f_p_lintri.m Given a mesh of triangles, this function computes several auxiliary variables and arrays (e.g., triangle areas, contribution of individual elements (triangles) to the stiffness/mass matrices and load vectors) to be used subsequently to assemble the various relevant stiffness and mass matrices as well as the load vectors used in the systems of ordinary differential equations. These arrays and variables are packed as fields into a structure “S” returned by the function. These arrays are formatted in a manner to optimize the generation (in terms of memory and speed) of the relevant matrices in sparse format (using the MATLAB function “sparse”). These matrices and load vectors are then computed from the structure “S” by the following functions. f_k_lintri_0.m Generates stiffness (conductance) matrix (K, in sparse format) from the structure “S” output by “f_p_lintri.m”, assuming that the stiffness parameter (e.g., extracellular conductance) is constant within every triangle. f_m_lintri_0.m Generates mass (capacitance) matrix (M, in sparse format) from the structure “S” output by “f_p_lintri.m”, assuming that the stiffness parameter (e.g., membrane capacitance) is constant within every triangle. f_f_lintri_0.m Generates load vector (F) from the structure “S” output by “f_p_lintri.m”, assuming that the load parameter (e.g., sodium current conductance per unit area) is constant within every triangle. f_k_lintri_1.m Generates stiffness (conductance) matrix (K, in sparse format) from the structure “S” output by “f_p_lintri.m”, assuming that the stiffness parameter (e.g., extracellular conductance) is a linear function within every triangle. f_k_lintri_1pt.m Variation of f_k_lintri_1.m. Generates stiffness (conductance) matrix (K, in sparse format) from the structure “S” output by “f_p_lintri.m”, assuming that the stiffness parameter (e.g., extracellular conductance) is a linear function within every triangle and there is another stiffness parameter that is uniform over each triangle. f_m_lintri_1.m Generates mass (capacitance) matrix (M, in sparse format) from the structure “S” output by “f_p_lintri.m”, assuming that the stiffness parameter (e.g., membrane capacitance) is a linear function within every triangle. f_f_lintri_1.m Generates load vector (F) from the structure “S” output by “f_p_lintri.m”, assuming that the load parameter (e.g., sodium current conductance per unit area) is a linear function within every triangle. Function used to dilate a meshed object f_Dilate.m Function used to identify the surrounding of objects in a two-dimensional finite element mesh of triangles, e.g., the peripheral rim of "dense extracellular proteins" with a "width" of N triangles (default value: 2 triangles) around Na+ channel clusters. Function used for modeling the Na+ current mhjLR1_livshitz.m This function computes and returns gating parameters for the three gates m, h and j of the sodium current according to Luo and Rudy (Circ Res 1991, 68:1501-1526) with modifications by Livshitz and Rudy (Biophys J 2009, 97:1265-1276). This function was already used in previous projects and uploaded previously on Zenodo (Ivanovic and Kucera, J Physiol 2021, 599: 4779-4811, DOI: 10.5281/zenodo.5226268; Cells 2022, 11, 3477, DOI: 10.5281/zenodo.7271839).