1"use strict";(self.webpackChunkdune_copasi_docusaurus=self.webpackChunkdune_copasi_docusaurus||[]).push([[9337],{1933:e=>{e.exports=JSON.parse('{"archive":{"blogPosts":[{"id":"cardiage-ep-v1","metadata":{"permalink":"/tutorials/cardiage-ep-v1","source":"@site/tutorials/cardiac_ep_1.mdx","title":"Cardiac Electrophysiology - I","description":"A tutorial for cardiac electrophysiology simulations using DuneCopasi.","date":"2024-03-19T20:55:29.000Z","tags":[{"inline":true,"label":"Electrophysiology","permalink":"/tutorials/tags/electrophysiology"},{"inline":true,"label":"Single Compartment","permalink":"/tutorials/tags/single-compartment"}],"readingTime":9.5,"hasTruncateMarker":true,"authors":[{"name":"Dylan Vermoortele","title":"Postdoctoral Researcher, Cardiovascular Imaging and Dynamics, KU Leuven, Belgium","url":"https://gbiomed.kuleuven.be/english/research/50000635/50508167/Team/00134802","page":{"permalink":"/tutorials/authors/dylan"},"socials":{"newsletter":"https://gbiomed.kuleuven.be/english/research/50000635/50508167/Team/00134802"},"imageURL":"https://gbiomed.kuleuven.be/english/research/50000635/50508167/Team/00134802/@@images/person_photo/thumb","key":"dylan"}],"frontMatter":{"title":"Cardiac Electrophysiology - I","description":"A tutorial for cardiac electrophysiology simulations using DuneCopasi.","slug":"cardiage-ep-v1","authors":["dylan"],"tags":["Electrophysiology","Single Compartment"],"hide_table_of_contents":false},"unlisted":false},"content":"import Tabs from \'@theme/Tabs\';\\nimport TabItem from \'@theme/TabItem\';\\nimport CodeBlock from \'@theme/CodeBlock\';\\nimport ExampleTerminal from \'@site/src/components/ExampleTerminal\'\\n\\nimport EPSimulation1INI from \'!!raw-loader!./assets/ini/cardiac_ep_1.ini\';\\n\\nIn this blog I will show you can configure the DuneCopasi simulation environment to perform a very simple cardiac electrophysiology simulation. I will show you how you can simulate the Mitchell-Schaeffer cardiomyocyte electrophysiology model in a 2D environment using the DuneCopasi Docker container.\\n\\n\\n\\n{/* truncate */}\\n\\nWe will simulate the Mitchell-Schaefer model where we stimulate the bottom left corner using a Gaussian point source. This will result in the activation and repolarization of our 2D square grid. Lets get started.\\n\\n## The mathematical model of cardiac electrophysiology\\n\\n### The Mitchell-Schaeffer model of cell electrophysiology\\n\\nCardiomyocytes maintain an ion concentration gradient between the intracellular and extracellular space. This ion concentration gradient results in a potential across the cell membrane (i.e. the transmembrane potential). A typical resting membrane potential is around $-85 mV$. A disturbance of this resting state can result in a cardiac action potential. This is a brief change in membrane potential of heart cells resulting in a typical transmembrane potential waveform. This typical waveform is determined by the dynamics of the ion channels in the cell membrane. In the very simple Mitchell-Schaeffer model the action potential is described by 2 ion channel currents. The governing equations are given by\\n\\n$$\\n\\\\begin{aligned}\\n C_m \\\\partial_t V_m = &J_{in}(V_m,h) + J_{out}(V_m),\\\\\\\\\\n &J_{in}(V_m,h) = \\\\frac{z {V_m}^2 (1-V_m)}{\\\\tau_{in}},\\\\\\\\\\n &J_{out}(V_m) = \\\\frac{V_m}{\\\\tau_{out}},\\n\\\\end{aligned}\\n$$\\n\\nwhere $\\\\tau_{in}$ and $\\\\tau_{out}$ are time-scale variables, $J_{in}$ and $J_{out}$ describe the respective inward and outward current, $V_m$ is the normalized transmembrane potential and $h$ is a gating variable. The dynamics of the gating variable are given by\\n\\n$$\\n\\\\frac{\\\\partial z}{\\\\partial t} =\\n\\\\begin{cases}\\n \\\\frac{z-1}{\\\\tau_{open}},& \\\\text{if } u\\\\leq u_0,\\\\\\\\\\n \\\\frac{z}{\\\\tau_{close}},& \\\\text{if } u > u_0.\\n\\\\end{cases}\\n$$\\nwhere $\\\\tau_{in}$ and $\\\\tau_{out}$ are time-scale variables for the gating variable and $u_0$ defines the critical threshold for depolarization.\\n\\n### The bi-domain equation of cardiac electrophysiology\\n\\nIn the section above we described a model capturing the basic electrophysiology of cardiomyocytes. Now we want to extend this by coupling multiple cardiomyocytes yielding a spatial description of the cardiac electrophysiology. To obtain this we will depart from Ohm\'s law stating\\n\\n$$\\n J = \\\\sigma E.\\n$$\\n\\nIn our specific case we have two compartments. We will use an intracellular and extracellular current density and potential. We obtain\\n\\n$$\\n\\\\begin{aligned}\\n J_{i} &= \\\\sigma_i E_i = \\\\sigma_i \\\\nabla \\\\varphi_i,\\\\\\\\\\n J_{e} &= \\\\sigma_e E_e = \\\\sigma_e \\\\nabla \\\\varphi_e. \\\\\\\\\\n\\\\end{aligned}\\n$$\\n\\nFurthermore, we assume that any current leaving the intracellular space flows into the extracellular (and vice-versa) and that there is no charge accumulation. Therefore, the change in current density in the extracellular space has to be equal to the current flowing into the intracellular space. Using Gauss\' law (i.e. the divergence theorem) this can be written as\\n\\n$$\\n \\\\nabla \\\\cdot J_{i} + \\\\nabla \\\\cdot J_{e} = 0 = \\\\nabla \\\\cdot \\\\sigma_i \\\\nabla \\\\varphi_i + \\\\nabla \\\\cdot \\\\sigma_e \\\\nabla \\\\varphi_e.\\n$$\\n\\nThe current flowing from the intracellular into the extracellular space is nothing else than the transmembrane current $I_m$ (per surface area). The conservation of current can be rewritten as\\n\\n$$\\n\\\\begin{aligned}\\n I_m &= \\\\nabla \\\\cdot \\\\sigma_i \\\\nabla \\\\varphi_i = -\\\\nabla \\\\cdot \\\\sigma_e \\\\nabla \\\\varphi_e.\\\\\\\\\\n\\\\end{aligned}\\n$$\\n\\nThe transmembrane current $I_m$ can be modeled using the cable equation where we combine the ion channel current $I_{ion}$ with a capacitative component\\n\\n$$\\n\\\\begin{aligned}\\n I_m &= A_m \\\\Big( C_m \\\\frac{\\\\partial V_m }{\\\\partial t} + I_{ion} \\\\Big)\\\\\\\\\\n\\\\end{aligned}\\n$$\\n\\nwhere $V_m = \\\\varphi_i - \\\\varphi_e$ is the transmembrane potential and $A_m$ is the surface density. Rewriting these equations (by eliminating $\\\\varphi_i$) yields the standard form of the bidomain equation of cardiac electrophysiology:\\n\\n$$\\n\\\\begin{aligned}\\n A_m \\\\Big( C_m \\\\frac{\\\\partial V_m }{\\\\partial t} + I_{ion} \\\\Big) &= \\\\nabla \\\\cdot \\\\sigma_i \\\\nabla V_m + \\\\nabla \\\\cdot \\\\sigma_i \\\\nabla \\\\varphi_e ,\\\\\\\\\\n 0 &= \\\\nabla \\\\cdot \\\\sigma_i \\\\nabla V_m + \\\\nabla \\\\cdot (\\\\sigma_i + \\\\sigma_e) \\\\nabla \\\\varphi_e . \\\\\\\\\\n\\\\end{aligned}\\n$$\\n\\n## EP simulations using DuneCopasi\\n\\nNow that we have discussed the mathematical problem we can tend our attention to implementing this so that we can simulate it using DuneCopasi. To do this we will write a [configuration file](/docs/param_tree) named `cardiac_ep_1.ini` that will entail all the information that is needed for the DuneCopasi solver to construct the problem at hand.\\nWe start by creating a working directory where we will save the simulation files.\\n\\n<ExampleTerminal prompt=\\"~/workdir$\\" title=\\"Creating a working directory\\">\\n {[\\n {\\n input: \\"mkdir -m o+rw workdir\\",\\n prompt: \\"~$\\",\\n },{\\n input: \\"cd workdir\\",\\n prompt: \\"~$\\",\\n }\\n ]}\\n</ExampleTerminal>\\n\\n### Setting up the Mitchell-Schaefer model\\n\\nLet us get started with creating the ini file for our electrophysiology simulation. Open your favorite text editor and
1you can start building the configuration file. First we will define our grid.\\n\\n#### The grid properties\\n\\nThe grid is defined as follows:\\n\\n\\n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - grid\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(0,4).join(\'\\\\n\')}\\n</CodeBlock>\\n\\n This will create a square grid existing of 2 triangles with 3 by 3 extension. Thereafter the grid will be refined 8 times. The solver can also load Gmsh meshes but we will come to that at a later point. For the moment this will be all that we need to define our grid. The next step is the ```parser_context```.\\n\\n#### The parser context\\n\\n In the parser context we will define the constants and functions that we will use in our model. This is easiest seen by looking at an example:\\n\\n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - parser context\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(5,36).join(\'\\\\n\')}\\n</CodeBlock>\\n\\n To create a constant variable lets say ```tau_in``` we simply write ```tau_in.type = constant```)). The value of the constant then assigned by writing ```tau_in.value = 0.1```. In a similar fashion a function is defined. For our Gaussian function which we name ```gauss``` we simply write ```gauss.type = function```. The function itself is then defined by typing ```gauss.expression = x,y,z : exp(...``` where the ```x,y,z``` before the colon denote the function arguments and the function itself is defined on the right of the colon. Note that standard mathematical functions such as ```exp```, ```sin```, etc. are available from the underlying parser allowing great freedom to define mathematical expressions. At this point we already have a grid and defined some constants and functions, now we will start defining the problem.\\n\\n#### Defining the compartments\\n\\nAlthough this example is fairly simple and we only need 1 compartment (which is the complete square), the concept to define multiple compartments is similar and therefore we will already discuss this briefly. Defining a compartment is similar to defining a defining a function. For our example we have one compartment which we will call the domain. To define this compartment type under the compartments section ```domain.type = expression```. This indicates that we will provide an expression which will be evaluated for each entity and if this expression evaluates to be non-zero then it will be added to the compartment. Here we only have one compartment which is all of the space so we set this expression equal to the constant 1.\\n\\n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - compartments\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(37,40).join(\'\\\\n\')}\\n</CodeBlock>\\n\\n#### The model properties\\n\\nThe first step in defining the model is to chose the parser that will be used to define the problem. Here we will go for [ExprTk](https://github.com/ArashPartow/exprtk).\\n\\n\\n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - model properties\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(41,43).join(\'\\\\n\')}\\n</CodeBlock>\\n\\nOnce that is selected we will start defining the first variable of our problem. In this case we will start with the transmembrane potential which we refer to as $$V_m$$. To do this we will create a subsection ```model.scalar_field.V``` with the properties of our problem.\\n\\n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - scalar field V\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(44,53).join(\'\\\\n\')}\\n</CodeBlock>\\n\\nThe ```compartment``` property will define the compartment on which the `scalar_field` $V$ is defined. Here we only had one compartment ```domain``` that corresponds to the complete grid. Since the equation determining $V$ is time-dependent we have to store the previous time-step for use in the time-stepping scheme. We do this by setting ```storage.expression``` to 1. The following 2 lines define the cross diffusion of respectively $V_m$ and $\\\\varphi_e$. Then we define the reaction term of the problem and the associated Jacobian entries of the reaction term. Note that we defined a function in the reaction term explicitly using the position $$(x,y,z)$$ and time $$t$$ to assemble the problem. At last we provide a function defining the initial state of the problem. In a similar fashion we define the extracellular potential field $\\\\varphi$.\\n\\n\\
1n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - scalar field phi\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(54,62).join(\'\\\\n\')}\\n</CodeBlock>\\n\\nThe $\\\\varphi$ equation is not time-dependent so we set ```storage.expression``` to zero. Further, we add the cross diffusion components and set the reaction term to zero. Note that explicitly setting the reaction term to zero is not necessary as the solver sets all non-defined properties to zero. Thereafter, we add a constraint to a small boundary region in the bottom left corner. Note that adding this constraint is technically not necessary, but since the $\\\\varphi_e$ variable only appears as gradients the solution is only determined up to a constant. So by setting the bottom left corner to zero we consider the potentials to be with respect to the potential in the bottom left corner. At last we provide the initial value. Similarly we define the gating variables $z$.\\n\\n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - scalar field z\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(63,69).join(\'\\\\n\')}\\n</CodeBlock>\\n\\nSince these equations has no diffusion components, we simply don\'t define them rendering the coefficients of the diffusion and cross-diffusion zero. By defining the problem as time-dependent with a reaction term and no diffusion this is in fact just a ordinary differential equation. So at this point we have defined the mathematical properties of the problem. The next step is to define the underlying time-step approach to solve these equations. This is done by setting the ```model.time_step_operator``` properties.\\n\\n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - time step operator\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(70,76).join(\'\\\\n\')}\\n</CodeBlock>\\n\\nFor this problem we will use an implicit-Euler time-stepping scheme and then we define the initial and maximal time-step. The simulation will step from time-point ```time_start``` to ```time_end```. In our case this will be from $0$ (default value) to $800$. To solve the linear problem we will use the BiCGSTAB iterative solver. The last step is to write the solution to a file. This is done by initializing the writer.\\n\\n<CodeBlock title=\\"~/workdir/cardiac_ep_1.ini - writer\\" language=\\"ini\\">\\n{EPSimulation1INI.split(/\\\\n/).slice(77,81).join(\'\\\\n\')}\\n</CodeBlock>\\n\\nAt this point we have defined all the properties of the problem and we are ready to simulate the model! Let us save the configuration file to ``cardiac_ep_1.ini``.\\n\\n<details>\\n<summary>\\nClick here to expand/collapse the complete configuration file\\n</summary>\\n<CodeBlock showLineNumbers title=\\"~/workdir/cardiac_ep_1.ini\\" language=\\"ini\\">\\n{EPSimulation1INI}\\n</CodeBlock>\\n</details>\\n\\n### Running the Mitchell-Schaefer model\\n\\nRunning the code is the same as described in the [Docker Documentation](/docs/install_docker#run-the-program), simply call:\\n\\n<ExampleTerminal prompt=\\"~/workdir$\\" title=\\"Running DuneCopasi (with docker)\\">\\n {[\\n {\\n input: \\"docker run -v $PWD:/dunecopasi dune-copasi --config=cardiac_ep_1.ini\\",\\n output: [\\n \\"[2024-03-21 12:38:32.479] [info] Reading configuration file \'~/workdir/cardiac_ep_1.ini\'\\",\\n \\"[2024-03-21 12:38:32.479] [info] Starting dune-copasi (version: 2.0.1)\\",\\n \\"[2024-03-21 12:38:32.580] [info] Axis names are set to: [\\\\\\"x\\\\\\", \\\\\\"y\\\\\\", \\\\\\"z\\\\\\"]\\",\\n \\"...\\",\\n \\"...\\",\\n \\"[2024-03-21 12:38:32.723] [info] dune-copasi successfully finished :)\\",\\n ]\\n }\\n ]}\\n</ExampleTerminal>\\n\\nThe simulation will start and write the ```*.vtk``` files to the specified directory. You can easily inspect these files by opening them in [Paraview](https://www.paraview.org/).\\n\\n\\n\\nAt this point we have successfully simulated the activation of a 2D square grid using the Mitchell-Schaefer cell model. Although this simulation model is very simple it does already illustrate the capabilities of the DuneCopasi solver. In the next post we will dig deeper into the available options and explore true multi-compartmental problems."}]}}')}}]);
Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.