**Simula SpringerBriefs on Computing** Reports on Computational Physiology **13**

**Kimberly J. McCabe** *Editor*

# **Computational Physiology** Simula Summer School 2022 − Student Reports

### **Simula SpringerBriefs on Computing**

### Reports on Computational Physiology

Volume 13

### **Series Editors**

Kimberly J. McCabe Simula Research Laboratory, Oslo, Norway Andrew D. McCulloch, Institute for Engineering in Medicine, University of California San Diego, La Jolla, California, USA; Department of Bioengineering, University of California San Diego, La Jolla, California, USA

### **Managing Editor**

Jennifer Hazen, Simula Research Laboratory, Oslo, Norway

### **Editor-in-Chief**

Aslak Tveito, Simula Research Laboratory, Oslo, Norway

### **About this series**

In 2016, Springer and Simula launched an Open Access series called the Simula SpringerBriefs on Computing. This series aims to provide concise introductions to the research areas in which Simula specializes: scientific computing, software engineering, communication systems, machine learning and cybersecurity. These books are written for graduate students, researchers, professionals and others who are keenly interested in the science of computing, and each volume presents a compact, state-of-the-art disciplinary overview and raises essential critical questions in the field.

Simula's focus on computational physiology has grown considerably over the last decade, with the development of multi-scale mathematical models of excitable tissues (brain and heart) that are becoming increasingly complex and accurate. This sub-series represents a new branch of the SimulaSpringer Briefs that is specifically focused on computational physiology. Each volume in this series will explore multiple physiological questions and the models developed to address them. Each of the questions will, in turn, be packaged into a short report format that provides a succinct summary of the findings and, whenever possible, the software used will be made publicly available.

Kimberly J. McCabe Editor

## Computational Physiology

Simula Summer School 2022 − Student Reports

*Editor*  Kimberly J. McCabe Simula Research Laboratory Oslo, Norway

ISSN 2512-1677 ISSN 2512-1685 (electronic) Simula SpringerBriefs on Computing ISBN 978-3-031-25373-7 ISBN 978-3-031-25374-4 (eBook) https://doi.org/10.1007/978-3-031-25374-4 ISSN 2730-7735 ISSN 2730-7743 (electronic) Reports on Computational Physiology

Mathematics Subject Classification (2020): 92-10, 92C99

© The Editor(s) (if applicable) and The Author(s) 2023. This book is an open access publication.

**Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

### **Preface**

Since 2014, we have organized an annual summer school in computational physiology. The school starts in June each year and the graduate students spend two weeks in Oslo learning the principles underlying mathematical models commonly used in studying the heart and the brain. At the end of their stay in Oslo, the students are assigned a research project to work on over the summer. In August the students travel to the University of California, San Diego to present their fndings. Each year, we have been duly impressed by the students' progress and we have often seen that the results contain the rudiments of a scientifc paper.

Starting in the 2021 edition of the summer school, we have taken the course one step further and aim to conclude every project with a scientifc report that passes rigorous peer review as a publication in this new series called *Simula SpringerBriefs on Computing – reports on computational physiology*.

One advantage of this course adjustment is that we have the opportunity to introduce students to scientifc writing. To ensure the students get the best introduction in the shortest amount of time, we have commissioned a professional introduction to science writing by Nature. The students participate in a 2-day *Nature Masterclasses* workshop, led by two editors from Nature journals, in order to strengthen skills in high quality scientifc writing and publishing. The workshop is tailored to publications in the feld of computational physiology and allows students to gather real-time feedback on their reports.

We would like to emphasise that the contributions in this series are brief reports based on the intensive research projects assigned during the summer school. Each report addresses a specifc problem of importance in physiology and presents a succinct summary of the fndings (8-15 pages). We do not require that results represent new scientifc results; rather, they can reproduce or supplement earlier computational studies or experimental fndings. The physiological question under consideration should be clearly formulated, the mathematical models should be defned in a manner readable by others at the same level of expertise, and the software used should, if possible, be made publicly available. All reports in this series are subjected to peer-review by the other students and supervisors in the program.

We would like to express our gratitude for the very fruitful collaboration with Springer -Nature and in particular with Dr. Martin Peters, the Executive Editor for Mathematics, Computational Science and Engineering.

The editors of *Simula SpringerBriefs on Computing – reports on computational physiology*:

Oslo, Norway December 2022

*Kimberly J McCabe Andrew D McCulloch Aslak Tveito Jennifer L Hazen*

### **Acknowledgements**

The Simula Summer School in Computational Physiology is a team efort, with many scientists contributing their time to give lectures and advise projects for the students. We would like to thank the lecturers and project advisors for their expertise and willing participation in the course: Dr. Hermenegild Arevalo, Dr. Gabriel Balaban, Dr. Jonas van den Brink, Dr. Francisco Contijoch, Dr. Cecile Daversin-Catty, Dr. Andrew ´ Edwards, Dr. Ada Johanne Ellingsrud, Dr. Nickolas Forsch, Dr. Ingeborg Gjerde, Dr. Ilse van Herck, Dr. Glenn Lines, Dr. Molly Maleckar, Dr. Kimberly McCabe, Dr. Andrew McCulloch, Lena Myklebust, Oscar Odeigah, Dr. Marie Rognes, Dr. Joakim Sundnes, Dr. Aslak Tveito, Julie Uv, Dr. Daniela Valdez-Jasso, Dr. Samuel Wall, and Kei Yamamoto.

Administrative support for the school was provided by Dr. Kimberly McCabe, Dr. Nickolas Forsch, Dr. Rachel Thomas, Dr. Stian Engen, and Hanie Tampus. At UCSD, we received vital administrative support from Jennifer Stowe and Marlene San Pedro. Formatting support for this volume was provided by Dr. Jennifer Hazen.

The Simula Summer School in Computational Physiology is supported through the Simula-UiO-UCSD Research PhD Training Programme (SUURPh), an endeavour funded by the Norwegian Ministry of Education and Research. Additional fnancial support is derived from SIMENT, an INTPART mobility grant from the Norwegian Research Council. The school also received funding from Digital Life Norway (DLN) to support student mobility.

### **Contents**

### **Chapter 1 An Automated Cardiac Constitutive Modelling Framework with Evolutionary Strain Energy Functions**

Kristin Ludwicki\*, Leto L Riebel\*, Sophia Ohnemus\*, Frida M E Westby, Nickolas Forsch, and Gabriel Balaban \* These authors contributed equally to this work

**Abstract** Heart disease is the leading cause of mortality worldwide. Many cardiac diseases are associated with altered elastic energy relations of the heart tissue. However, the strain energy functions describing these characteristics are limited since they were designed manually for highly specifc experimental setups. In this study, we develop CHESRA (Cardiac Hyperelastic Evolutionary Symbolic Regression Algorithm), an automated constitutive modelling framework, to derive cardiac elastic strain energy functions directly from experimental data. Our results indicate that CHESRA fnds functions that reproduce mechanical tissue properties from experimental data whilst controlling function complexity. Our novel approach has the potential to fnd strain-energy functions that ft to various experimental data sets and may contribute to automatically building mathematical models to understand clinical observations of heart diseases.

Sophia Ohnemus Institute for Experimental Cardiovascular Medicine, University Heart Center Freiburg - Bad Krozingen, Medical Center and Faculty of Medicine, University of Freiburg, DE

Frida M E Westby Department of Informatics, University of Oslo, NO, Department of Computational Physiology, Simula Research Laboratory, NO

Nickolas Forsch Department of Computational Physiology, Simula Research Laboratory, NO

Gabriel Balaban (corresponding author) Department of Computational Physiology, Simula Research Laboratory, NO e-mail: gabrib@ simula.no

Kristin Ludwicki

Department of Physiology and Biophysics,Weill Cornell Medicine, New York, USA

Leto L Riebel Department of Computer Science, University of Oxford, UK

### **1.1 Introduction**

Cardiovascular disease is the leading cause of mortality worldwide, accounting for around 19 millions deaths in 2020 [1]. In many cardiac diseases, the mechanical and structural properties of the heart are altered. Examples include hypertrophic cardiomyopathy, which is characterized by a thickening of the heart walls [2] and dilated cardiomyopathy, a disease which causes the ventricular walls to thin and stretch [3]. Another example is myocardial infarction, commonly known as a heart attack, during which oxygen and nutrient supply to the cardiac tissue is interrupted [4]. Myocardial infarction results in altered mechanical loading conditions [5], structural changes, increased tissue stifness, and ultimately ventricular dysfunction [6]. A better understanding of these pathological mechanisms may improve clinical diagnosis and treatment options.

To characterize the elastic properties of the heart, stress-strain relations of the myocardium have been extensively studied experimentally. Assessment of passive shear properties of pig ventricular myocardium by Dokos et al. [7] indicates that the myocardium presents as an orthotropic material since the shear response difers along three mutually orthogonal directions. Demer and Yin [8] and Yin et al. [9] measured the stress-strain relationship of passive non-contracting canine myocardium for simultaneous biaxial stretching and observed highly nonlinear, anisotropic, and viscoelastic behavior. Novak et al. [10] performed more biaxial experiments in canine myocardium from diferent regions of the left ventricle and found that qualitatively the mechanical properties were similar, but quantitatively they difered. For example, samples from sub-endocardium and sub-epicardium of the left ventricular free wall tended to be more rigid than samples from the mid-myocardium.

In order to mathematically represent the mechanical properties of the heart tissue, a number of strain-energy functions (SEF) have been proposed. These SEF describe the potential energy density of the heart tissue depending on its deformation and have been developed based on a few key considerations. Considerations include structural metrics and other experimental results, as well as assumptions on physical properties of the system. Although the myocardium appears to be a viscoelastic material, the relaxation time of the viscoelastic response is long compared to the cardiac cycle [11]. Therefore, most SEF, with exceptions as presented in [12], model it as hyperelastic for simplicity reasons. Examples for SEF based on the assumption of transverse isotropy are given in [13, 14], whereas [11, 15, 16, 17] introduce orthotropic models. However, creating these SEF is labor intensive as they must be derived manually and, thus, only a limited amount of experimental data can be considered.

Evolutionary algorithms have the potential to automatically fnd functions that describe experimental data with minimal human guidance. This has been applied already to various research felds; for example, Doglioni et al. used evolutionary polynomial regression (EPR) to build a function for the relationship between air and water temperature [18]. Their method not only considered suitability of the discovered function to describe the given data, but also aimed to minimise the complexity of the formula. Montes et al. used a modifed EPR framework to create self-cleansing models for new sewer systems [19] and several diferent machine learning methods, including EPR, were used by Jamei et al. to model the density of hybrid nanofuids [20]. Furthermore, Javadi et al. established an EPR method to build constitutive equations describing the deformation of bridges and tunnels [21].

In this work we present the development of a Cardiac Hyperelastic Evolutionary Symbolic Regression Algorithm (CHESRA) that derives hyperelastic SEF for ventricular myocardium. In contrast to EPR algorithms, CHESRA is not restricted to polynomials, but searches the whole space of mathematical functions for SEF describing the given experimental data. Our results highlight the potential to automatically fnd SEF which reproduce experimental fndings of cardiac elastic properties of diferent species.

### **1.2 Methods**

In the following subsections we outline the individual components of CHESRA. First, in 1.2.1 we briefy introduce the physical principles of hyperelastic continuum mechanics that we use to defne the SEF. Then, in 1.2.2 we outline how we represent the SEF computationally as function trees. Subsequently, we explain how we ft the SEF to experimental data in 1.2.3. Lastly, we go through the individual steps of the evolutionary symbolic regression (ESR) algorithm in detail in 1.2.4.

### **1.2.1 Basic Principles of Hyperelastic Continuum Mechanics**

In order to defne our SEF, we build upon the framework developed by Holzapfel and Ogden for cardiac hyperelasticity [11]. Here, the structure of the myocardium is characterized by three orthonormal basis vectors: the fbre axis associated with the prevailing cell orientation, the sheet axis defned as the direction in the plane of the cell sheets perpendicular to the fbre direction, and the sheet-normal axis perpendicular to and .

Any point within the material can be described by a vector **X** = Í = ,, **e** , where **e** denote the basis vectors of the , , and axis. Upon deformation, the position of each point changes and the new position can be described by a vector **x**. The fundamental quantity describing such a deformation is the deformation gradient **F**, defned by = / with , = , , and where = det**F** = 1 for an incompressible material such as the myocardium. The right Cauchy-Green tensor

$$\mathbf{C} = \mathbf{F}^T \mathbf{F},\tag{1.1}$$

and the Green-Lagrange strain tensor

$$\mathbf{E} = \frac{1}{2}(\mathbf{C} - \mathbf{I}),\tag{1.2}$$

are associated with **F**, with **I** being the identity tensor.

Since these quantities do not only account for deformation but also for rotation and translation, in the following, we consider the invariants of **C** which are defned as

$$\begin{aligned} I\_1 &= \text{tr}\mathbf{C}, & I\_2 &= \frac{1}{2} \left[ I\_1^2 - \text{tr}\mathbf{C} \right] & & I\_3 &= \text{det}\mathbf{C}, \\ I\_{4i} &= \mathbf{e}\_i^T (\mathbf{C}\mathbf{e}\_i), & I\_{5i} &= \mathbf{e}\_i^T (\mathbf{C}^2\mathbf{e}\_i), & & I\_{8if} = \mathbf{e}\_i^T (\mathbf{C}\mathbf{e}\_j), \end{aligned} \tag{1.3}$$

with ≠ ∈ { , , }. Defning a SEF based on these invariants ensures objectivity. The SEF proposed by Holzapfel and Ogden [11] is given by

$$\psi = \frac{a}{2b} \exp\{b(I\_1 - 3)\} + \sum\_{i=f,s} \frac{a\_i}{2b\_i} \left\{ \exp\{b\_i(I\_{4i} - 1)^2\} - 1 \right\}$$

$$+ \frac{a\_{fs}}{2b\_{fs}} \left\{ \exp\{b\_{fs}I\_{8fs}^2\} - 1 \right\}. \tag{1.4}$$

Here, , , , , , , , and are material parameters chosen to either ft the corresponding Cauchy stress tensor

$$\boldsymbol{\sigma} = \boldsymbol{J}^{-1} \mathbf{F} \sum\_{k} \frac{\partial \boldsymbol{\psi}}{\partial I\_{k}} \frac{\partial I\_{k}}{\partial \mathbf{F}} \tag{1.5}$$

to the shear data from Dokos et al. [7], or the associated second Piola-Kirchhof stress tensor

$$\mathbf{S} = J \mathbf{F}^{-1} \sigma \mathbf{F}^{-T} \tag{1.6}$$

to the biaxial stretch data by Yin et al. [9] (see Section <IP_ADDRESS> for details about these experimental data sets).

### **1.2.2 Representation of the Strain-Energy Functions as Function Trees**

In CHESRA, SEF are represented as a list of nodes composing function trees. Details on function tree implementation and use cases can be drawn from [22, 23]. In our implementation, each node in the function tree has a type, a value, up to one parent, and up to two children (one left and one right). Only the root node of the function tree has no parent, further explained in <IP_ADDRESS>. Nodes with no children are called 'leaves'. Sub-trees within the function tree represent individual terms of the SEF.

Nodes are either a symbol, an operand, or a pre-operand. Nodes of type symbol can either be a constant material parameter or an invariant. Material parameters are represented by the symbols {1, 2, 3, 4, 5, 6, 7, 8}. These material parameters are later replaced by a constant real number to ft the SEF to experimental data, as described in 1.2.3. Similarly, invariant nodes are assigned a symbol from the set of invariants {1, 2, 3, 4 , 4, 4, 5 , 5, 5, 8 , 8 , 8 }, as defned in Equation 1.3. Operand nodes can either have the value + (addition), − (subtraction), ÷ (division) or ∗ (multiplication); and pre-operands are either (to the power of 2), (euler's number '' to the power of) or - (negation). While operands require two terms, for example + , a pre-operand requires only one, for instance , which is written before the term in many coding languages, such as in python: (,2). An example function tree is shown in Figure 1.1.

Fig. 1.1: Example SEF and corresponding function tree with node types indicated. The function tree is converted to a real function by applying an in-order traversal.

If a node of type symbol, *i.e.* either a material parameter or an invariant, has a parent node, it has to be either an operand or a pre-operand to be a valid function. Similarly, symbol nodes can only have children of type operand. Pre-operands are applied to terms, and hence can only have exactly one child, which has to be of type symbol or pre-operand. Their parent has to be either an operand or another pre-operand. Operands have a parent and exactly one child, both of which cannot be another operand. If an operand is a right child, its child needs to also be a right child and vice versa, to avoid two symbols next to each other. These rules are summarised in Table 1.1. We apply brackets only around the terms pre-operators are applied to, all other operations are carried out according to their priority, e.g. ∗ before +.

Using the rules summarized in Table 1.1, to convert our function trees into functions we apply an in-order traversal, where any node is read out as: value of its left child, its own value, value of its right child. Figure 1.1 shows an example of a function and one of its possible corresponding function trees, depending on the order the terms are created in. Note that since our algorithm does not perform any searches within the tree, there is no need for it to be balanced.

Representing SEF as function trees has several benefts. Firstly, individual terms can be identifed as sub-trees, allowing for pre-operands to be applied to parts of the equation only. This would be more difcult, for example, in a linked list structure. Secondly, terms can be manipulated, swapped, and deleted easily during mutation and mating, which are explained in <IP_ADDRESS>.

Table 1.1: Summary of node types, possible node values, and node connection rules within the function tree.

### **1.2.3 Fitting the Strain-Energy Functions to Experimental Data**

For each SEF we ft the material parameters to experimental data sets in a separate step, in order to ensure that only the material parameters difer for diferent experimental setups while the general form of the SEF is the same. In this work we consider the shear data of Dokos et al. [7] and the biaxial stretch data by Yin et al. [9], which we briefy describe in the following.

#### **<IP_ADDRESS> Experimental Data under Consideration**

Dokos et al. [7] measured the shear stress, corresponding to the Cauchy stress tensor components versus the amount of shear in a cube of pig left ventricular myocardium sheared in the , , or plane (see Figure 1.2a). They considered diferent shear modes ( ) referring to shear in the plane in direction with ≠ ∈ { , , }.

The biaxial stretch data from Yin et al. [9] was collected in canine left ventricular tissue which was stretched simultaneously in the fbre and sheets direction. Figure 1.2b shows the measured second Piola-Kirchhof stress versus strain and respectively versus for three diferent constant strain ratios = /.

For both data sets, we used the digitizer software WebPlotDigitizer [24] to digitize the data from the graphs shown in [7] and [9].

#### **<IP_ADDRESS> Fitting to the Experimental Data**

In order to ft the material parameters of a SEF to the shear data set we calculated the predicted Cauchy stress tensor components according to Equation 1.5. Then, we used the Python package lmft [25] to minimize the residuals,

$$res\_k = s\_{if}(\gamma\_k) - \sigma\_{if}(\gamma\_k),\tag{1.7}$$

Fig. 1.2: Experimental data sets under consideration. (a) Plot of the digitized data collected by Dokos et al. [7] in pig left ventricular myocardium for shear stress versus the amount of shear for shear ( ) in the plane in direction with ≠ ∈ { , , }. Results for ( ) and () shear are identical. (b) Visualization of the data by Yin et al. [9] for biaxial loading in the plane of canine left ventricular myocardium. The left plot shows the stress versus strain in the fbre direction and the right plot versus in the sheets direction for three diferent strain ratios = / = 2.05 (triangles), 1.02 (squares) and 0.48 (circles). In both cases we used WebPlotDigitizer [24] to digitize the data shown in [7] and [9].

between experimental data points and the corresponding calculated stress tensor component for all shear modes ( ) and all amounts of shear in the data set simultaneously.

Analogously, for ftting the material parameters to the biaxial data set we calculated the diagonal components of the second Piola-Kirchhof stress tensor according to Equation 1.6. We used the same Python package for simultaneously minimizing the residuals

$$res\_k = s\_{i\bar{i}}(E\_{i\bar{i},k}) - S\_{i\bar{i}}(E\_{i\bar{i},k}),\tag{1.8}$$

with = , , for all Green-Lagrange strain tensor components , and strain ratios = / in the data set.

### **1.2.4 Evolutionary Symbolic Regression Algorithm**

To develop CHESRA, we build upon the traditional ESR framework. First, populations of random SEF are created (see Section <IP_ADDRESS>). Then, every equation within each iteration, or generation, is assigned an error-score based on a user-designed ftness function. Our ftness function assesses the complexity of the SEF and how well it fts to the shear and/or biaxial data set (see Section <IP_ADDRESS>). The individuals with the best ftness are then selected, mated to generate a new population, and mutated to ensure diversity for the next generation (see Section <IP_ADDRESS>). This multi-step process continues until the maximum number of generations is reached. A schematic representation of CHESRA is illustrated in Figure 1.3.

Fig. 1.3: The main workfow of CHESRA.

### **<IP_ADDRESS> Setup and Initialisation**

To initialise SEF and their corresponding function trees, a maximum function length, set of invariants, and material parameters are defned. First, our algorithm creates a root node, which is always an invariant or a material parameter. Second, it picks a random length within the given maximum and extends the function tree until this chosen length is reached. The function can either be extended by adding a pre-operand in front of a term, or by attaching an operand and a new symbol (i.e., an invariant or material parameter) to an existing symbol. For the latter, the selected existing symbol may have at most one child. Details concerning function development can be found in Section 1.2.2.

#### **<IP_ADDRESS> Mutation and Mating of Strain-Energy Functions**

We mutate SEF by randomly selecting a node from the corresponding function tree to be mutated. The selected node is then replaced by a randomly chosen element of the same type (see 1.4).

To mate two parent SEF and create two new children, our algorithm creates a deep copy of each parent function tree, chooses one random sub-tree per copy, and then swaps them over. To make sure the resulting children SEF are valid functions, we restrict the chosen sub-trees to be proper terms and hence start with a symbol or pre-operand. Alternatively, one may also choose to ensure the sub-trees to be swapped to start with a node of the same type. Our chosen sub-trees may be of diferent size and hence children SEF may be of diferent lengths than their parents (see 1.5).

Mutation and mating are important to increase population diversity in CHESRA. However, too much diversity could also prohibit convergence towards an optimal solution. Therefore, we set rates that will dictate the probability of a given individual to mate, mate, and/or mutate, mutate. These rates were set to mate = 20% and mutate = 80% in all CHESRA experiments.

Fig. 1.4: An examplary function tree representing a SEF is mutated by choosing a random item and replacing it with a random value of the same type, as highlighted by the orange boxes.

#### **<IP_ADDRESS> Scoring Fitness of Strain-Energy Functions**

We designed a ftness function to ensure that CHESRA can penalize against poor ftting SEF with high complexity. Therefore, each SEF, or individual, in each generation is evaluated based on its complexity and how well it reproduces the results of a specifc experimental data set. It is this evaluation that allows the best individuals to be selected for mating, mutation, and population of the following generation.

Fig. 1.5: Mating between two parent SEF to generate two children SEF is implemented by swapping two random subtrees from the corresponding function trees as highlighted in the orange boxes.

To evaluate the ftness of each equation, frst we quantify its complexity as the sum of the function length eq (i.e., the number of nodes it consists of), and the number of pre-operators pre-op that it contains. The pre-operator count allows for penalization against nestedness, or the creation of functions with pre-operators within themselves.

$$Complexity = l\_{\text{eq}} + n\_{\text{pre-op}} \,. \tag{1.9}$$

Second, we assess how well each individual SEF reproduces experimental data. For this, the material parameters of each equation are frst ftted to the experimental data as described in <IP_ADDRESS>. Then, for each equation the sum of squared errors (SSE) is calculated as the sum of the squared residuals between the experimental data points and the ft,

$$SSE = \sum\_{k} res\_k^2,\tag{1.10}$$

with as defned in Equation 1.7 for the shear data or in Equation 1.8 for the biaxial stretch data.

The fnal ftness score is then given by

$$Fitness = (\alpha \* Compplexity) + SSE,\tag{1.11}$$

where the hyperparameter ensures appropriate balance between ftting and complexity. Numerical experiments conducted to fnd a suitable can be found in Section 1.3.1.

Where CHESRA is given data sets, the ftness function can be extended to

$$Fitness = (\alpha \* Compplexity) + \sum\_{N} SSE\_N. \tag{1.12}$$

Hence, in the following section 1.3.1 we used Equation 1.11 in order to fnd a SEF that can reproduce the shear data alone, while in section 1.3.2 we used Equation 1.12 to consider both the shear and the biaxial data set.

### **1.3 Results**

### **1.3.1 Trial 1: SEF Describing the Shear Data**

In our frst trial, we ran CHESRA while ftting the SEF to the shear data set from Dokos et al. [7]. We chose the hyperparameter of the ftness function (see Equation 1.10) by running CHESRA with = 10−5 , 10−3 , 10−1 , 1, 2, and 5 for 50 generations and 100 individuals each. 1.6 indicates that a value of = 0.1 lead to the best result as it allowed for SEF with low SSE and minimal complexity. Using this hyperparameter, the ftness score of the best SEF gradually decreased over the generations until a plateau was reached (see 1.7).

In this run, the SEF with the lowest ftness score derived using CHESRA reproduces the shear data well (see 1.8a) and is given by

$$\psi\_1 = \exp\left(p\_1 I\_{\mathbb{S}s} + p\_2 I\_{\mathbb{S}f} + 2I\_2 - I\_{4n} - p\_3\right). \tag{1.13}$$

In order to assess whether CHESRA is reproducible, we repeated this trial twice, with the best SEF given by

$$
\psi\_2 = -\left(I\_{\mathbb{S}f} - p\_1\right)^2 + p\_2 I\_{\mathbb{S}f} \left(\left(I\_{\mathbb{S}} + p\_3\right)^2 + I\_{4s}\right)^2,\tag{1.14}
$$

$$\psi\_3 = \exp\left(p\_1^2 I\_{\xi f} + (p\_2 + p\_3 - 1)I\_{8fs} + p\_4 I\_{4f} + I\_1 + 2I\_2 - I\_{4n} + p\_S\right). \tag{1.15}$$

The equations derived from the three trial runs have unique material parameters that can be found in 1.2. 1.8 shows a good ft to the shear data for all three SEF found.

Table 1.2: Material parameters to ft the SEF defned in Equation 1.13 (trial run 1), 1.14 (trial run 2), 1.15 (trial run 3) to the shear data set by Dokos et al. [7]

Fig. 1.6: SSE vs. complexity for diferent hyperparameters in the ftness function (see Equation 1.11). Six runs of CHESRA were completed in which each was run with one of the unique values listed in the legend. For experiment 1, we chose = 0.1 since it lead to SEF with low SSE and function complexity.

Fig. 1.7: Evolution of 100 individual SEF for 50 generations. Each point represents the ftness of the best SEF from each generation.

Fig. 1.8: In our numerical experiment 1, CHESRA found a SEF that reproduces the shear data of Dokos et al. [7] in all 3 trial runs. Markers represent the experimental data points whereas solid lines show the ft of the SEF derived using CHESRA.

### **1.3.2 Trial 2: SEF Describing the Shear and Biaxial Data**

As a next step, we extended CHESRA to fnd a SEF that reproduces the results of multiple experimental data sets. Therefore, we ran CHESRA with the same setup as in trial 1, but updated the ftness function to incorporate the SSE of both the shear and the biaxial data set, as described in Equation 1.12. Moreover, in this numerical experiment we used = 0 in order to assess quality of ft without penalizing for equation complexity. In addition, = 0 allowed convergence in a shorter number of generations as ftting to two diferent data sets is a challenging optimization problem.

The SEF found using CHESRA is given by

$$
\psi = \left( p\_1 p\_2 I\_{4f} - I\_{4f}^2 + \frac{p\_1}{I\_{5n}} - I\_{4s} \left[ p\_1 p\_2 I\_{5n} I\_{4f} - p\_3 p\_4^2 I\_1 + p\_5 + p\_6 + p\_7 \right]^2 \right)
$$

$$
$$

$$
\left( \begin{array}{c} \\ \end{array} \right)^2, \tag{1.16}
$$

where the material parameter values are provided in Table 1.3. Figure 1.9 shows that the equation derived in experiment 2 has a better ft to the biaxial data than to the shear data.

Table 1.3: Material parameters for ftting the SEF defned in Equation 1.16 to the shear and biaxial data sets.

Fig. 1.9: In experiment 2, CHESRA was extended to fnd a single SEF aimed to reproduce both the shear data from Dokos et al. (a) and biaxial stretch data from Yin et al. (b). Solid lines show the ft of the SEF found using CHESRA to the experimental data represented by the markers.

### **1.4 Discussion**

CHESRA automatically generates functions which reproduce experimental data of cardiac hyperelastic properties. The results of trial 1 indicate that CHESRA can derive SEF that ft the shear data from Dokos et al. [7] well. However, the quality of ft was not equal for all shear directions. Specifcally, in all trial runs the derived SEF did not ft to the shear data for the (nf) and (ns) direction, as shown in Figure 1.8. As the magnitude of shear stress is much lower in those cases, it is likely that the ftness function neglects these shear modes. Previous studies have used coefcient of determination in the ftness function [18, 21], which will be a key consideration in our future work.

To extend CHESRA to multiple experimental data sets, we included both shear [7] and biaxial data [9] in trial 2. Here, the biaxial data was ftted at the expense of the shear data. A possible cause for this could be a non-optimal individual and generation number. Since more data sets were given to the algorithm, optimization is challenging in only 50 generations. Therefore, increasing the number of generations and individuals may allow for better convergence and a more universal SEF to be found.

Ensuring that CHESRA derives simple equations was important to this study. Therefore, we included the hyperparameter to the ftness function to ensure penalization against long equations with nesting. While this is a valuable proof-of-concept, it has limitations. Specifcally, it does not take computational time into consideration. Long but simple equations might be penalized more compared to shorter, but more complex or nested functions. The time it takes for an equation to be calculated may be a better measure of simplicity and will be explored in next steps.

Moreover, the current version of CHESRA can only generate SEF using quadratic and exponential functions; however, this could be widened to a larger function space. In addition, CHESRA may be extended to even more experimental data sets. Future work may also investigate to what extent the functions generated by CHESRA fulfll the mathematical requirements of SEF, such as convexity and strong ellipticity [11]. Lastly, it would be interesting to compare our method to a Multivariate Adaptive Regression Spline (MARS) model since studies by Jamei et al. [20] suggest that this can be more accurate than ESR.

### **1.5 Conclusions**

In this work we present CHESRA, an ESR framework for automatically deriving strain-energy functions of the myocardium. Our algorithm allows for controlling function complexity and inputting multiple experimental data sources.

### **References**

regression–multi-objective genetic algorithm strategy. *Urban Water Journal*, 17(2):154–162, 2020.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 2 Electromechanical** *In Silico* **Testing Alters Predicted Drug-Induced Risk to Develop Torsade de Pointes**

Anna Busatto, Jonathan Krauß, Evianne Kruithof, Hermenegild Arevalo, and Ilse van Herck

**Abstract** Torsade de Pointes (TdP) is a type of ventricular tachycardia that can occur as a side efect of several medications. The Comprehensive *in vitro* Proarrhythmia Assay (CiPA) is a novel testing paradigm that utilizes single cell electrophysiological simulations to predict TdP risk for drugs that could potentially be used clinically. However, the efects on mechanical performance and mechano-electrical feedback are neglected. Here, we demonstrate that including electromechanical simulations in CiPA testing can provide additional insights into the predicted drug-induced TdP risk. In this work, we analyzed six drugs, namely fecainide, ibutilide, metronidazole, mexiletine, quinidine and ranolazine. We compared previously classifed risks (low, intermediate, high) with our fully coupled electromechanical simulation results based upon the action potential, the electromechanical window, and the maximum active tension [1]. For ranolazine and metronidazole the predicted risk changed from low to intermediate and intermediate to high, respectively. For the latter, while electrophysiological markers indicated a low risk, the active tension decreased by 58% which can result in a loss of heart function. Therefore, adding mechanics to CiPA testing results in an altered prediction of drug-related TdP risk.

Anna Busatto

Jonathan Krauß Institute of Biomedical Engineering, Karlsruhe Institute of Technology, DE

Evianne Kruithof Eindhoven University of Technology, NL

Hermenegild Arevalo Department of Computational Physiology, Simula Research Laboratory, NO

Ilse van Herck (corresponding author)

Department of Biomedical Engineering, University of Utah, USA, Scientifc Computing and Imaging Institute, University of Utah, USA

Department of Computational Physiology, Simula Research Laboratory, NO e-mail: ilse@ simula.no

### **2.1 Introduction**

Torsade de Pointes (TdP) is a form of abnormal heart rhythm often preceded or caused by a prolonged QT-interval of the ECG [2]. Due to the high risk of unwanted and dangerous side-efects, promising new drugs with a potential TdP risk have been excluded from the research and drug-development pipeline. This, in combination with the long and complicated process of drug approval, led to a lack of new drugs on the market for cardiac disease. Therefore, the Comprehensive *in vitro* Proarrhythmia Assay (CiPA) is performed to better predict the TdP risk for drugs that could potentially be used in the clinic. However, some of the tests can result in inaccurate predictions of the actual outcomes and lead to restrictive results, limiting the number of drugs allowed on the market.

There has been ongoing research on improving these models to increase the predictability of drug-induced efects in the heart [3]. Llopis *et al.* proposed CiPA testing using a population of models approach to obtain more accurate risk predictions by simulating the efects of altered ion channel functions on the electrophysiological behavior of cardiac cells [1].

The electromechanical window () is a biomarker for TdP risk introduced by Passini *et al.* for *in silico* drug testing [4]. Clinically, is defned as the diference between the duration of electrical and mechanical systole [4]. However, for the *in silico* test, they defned the as the time diference between calcium transient duration at 90% repolarization ( 90) and action potential duration at 90% repolarization (90) [4]. This assumption makes the simulation computationally less expensive and complex, but this simplifcation reduces accuracy in the mechanical component.

In previous work, the efects of electromechanical feedback are neglected and no insight into the mechanical function is obtained. This feedback alters the electrophysiological results, especially the calcium transient, and afects the biomarkers used to predict TdP risk. However, the efect of drugs on the mechanical performance can be investigated with electromechanical simulations. Therefore, in this work, a fully coupled electromechanical simulator SimCardEMS is used to perform CiPA tests [5]. Here, the cardiac mechanical function can be analyzed in terms of active tension, tissue deformations and where the latter can be altered due to electromechanical feedback. We performed several tests to compare our predicted TdP risks to previous classifcations such as the ones found in Llopis *et al.* [1].

### **2.2 Methods**

To assess TdP risk, we focused on selected biomarkers, namely maximum active tension a,max [kPa], 90 [ms], 90 [ms], and the [ms]. Using these metrics, we analyzed six diferent drugs from CiPA: fecainide, ibutilide, metronidazole, mexiletine, quinidine, and ranolazine with known cardiac efects. Each of the drugs has been classifed as high, intermediate, or low risk for TdP in Llopis *et* *al.* [1]. Here, the indicated risks are compared with the risks predicted by a fully coupled electromechanical model to investigate whether the mechanical biomarkers should also be taken into account for risk classifcation.

### **2.2.1 Drug Implementation**

The cardiac efects of the drugs are tested for three diferent plasma concentrations (). For each of the concentrations, the ion channel conductances () are multiplied with scaling factors () representing the efect of the tested drug. The values are calculated using the drug , the inhibitory concentration (50) and Hill coefcient (ℎ) value of the targeted ion channel:

$$SF = \frac{\text{g}\_{\text{drug}}}{\text{g}} = \left[1 + (\frac{PC}{IC\_{S0}})^h\right]^{-1} \tag{2.1}$$

Table 2.1: Ion channel scaling factors (SF) for quinidine using three diferent plasma concentrations (). Scaling factors for seven ion channels were calculated with Equation 2.1.

Equation 2.1 is derived from the standard Hill equation [1]. The , 50 and ℎ values are taken from Llopis *et al.*. The values are calculated for the considered ion channels Kr, Na, NaL, CaL, Ks, K1 and to at three diferent s which are half (0.5), normal (1.0) and double (2.0) the efective free therapeutic plasma concentration [1]. As an example, specifc values used to simulate the efect of quinidine are shown in Table 2.1.

### **2.2.2 SimCardEMS Simulations**

The calculated values are fnally used as input for the SimCardEMS solver, which is implemented in FEniCS [6], to simulate the efects of drugs on an endocardial tissue block. We solved the electrophysiological O'Hara-Rudy model which is strongly coupled with the mechanical Land model via intracellular calcium concentration and stretch activated channels [7, 8]. The active tension was scaled by a factor of fve. The Holzapfel model was used to describe the tissue material properties, where we increased the stifness by a factor of ten [9].

The simulations were performed on a 6x3x3 mm mesh representing an endocardial tissue slab with 1 mm and 0.5 mm spatial resolution for the mechanical and electrophysiological parts, respectively (Fig. 2.1). One corner of the mesh was fxed, and each of the connected planes were fxed in their respective normal directions, allowing the simulation of free contraction. The fbers were aligned in the longitudinal direction of the mesh and the tissue was activated by stimulating the entire mesh.

To reach steady state, 40 beats of 1000 ms each were simulated and used as control. Thereafter, the drugs were added into the simulation by loading the calculated s. Each of these simulations was run either for fve or ten beats (1000 ms each) depending on when steady state was reached; the goal was to investigate the efect of each tested drug on the electrophysiological and mechanical behavior of the cells. The biomarkers used for the analysis were extracted from the center of the mesh for the last beat of each simulation.

### **2.2.3 Evaluation**

The TdP risk classifcations were based upon maximum active tension a,max, 90, 90, and where

$$EM\_{\mathcal{W}} = CaTD\_{\mathcal{W}0} - APD\_{\mathcal{W}0} \tag{2.2}$$

Fig. 2.1: (Left) Mechanical mesh 6x3x3 mm with 1 mm resolution. (Right) Electrophysiological mesh 6x3x3 mm with 0.5 mm resolution. Three planes are fxed in their normal direction to allow free contraction of the material.

The metrics used for evaluation are visualized in Fig. 2.2. TdP risk was classifed as high or intermediate respectively in the following situations:

If none of these criteria were met, the drug was considered safe and classifed as low risk.

Fig. 2.2: Transmembrane voltage, calcium concentration, and active tension for control. The metrics used to evaluate the results are indicated in the fgure with dashed lines.

### **2.3 Results**

The biomarkers extracted for all simulations are shown in Table 2.2. For the control simulation we found a , of 55.0 kPa, 90 of 258 ms and of 134 ms. The associated risk classifcation thresholds for , were 44.0 kPa (intermediate) and 27.5 kPa (high). For 90, the thresholds were 273 ms (intermediate) and 284 ms (high), while for they were 121 ms (intermediate) and 107 ms (high).

### **2.3.1 Action Potential Duration**

The resulting action potentials for control as well as the six examined drugs at 1.0 are shown in Fig. 2.3. Ibutilide, quinidine and fecainide have 90 values of 597 ms, 446 ms and 352 ms respectively which are above the threshold for high risk drugs. Ranolazine has an 90 of 283 ms at 1.0 which fts in the range of intermediate risk drugs. 90 values for metronidazole and mexiletine are below

Table 2.2: Overview of the extracted biomarkers, a,max [kPa], 90 [ms], 90 [ms], [ms] and TdP risk classifcations for all performed simulations. 0.5, 1.0 and 2.0 results are given in this order for each of the six tested drugs.

the intermediate threshold and should therefore be classifed as low risk if only the 90 is considered.

All the chosen drugs were then examined for 0.5 and 2.0 as well. An example of the diferences of action potential shape for quinidine based on the PC is shown in Fig. 2.4. A concentration dependent efect is observed for 90 as it is the smallest for 0.5 and increases with an increase in PC. An overview of the 90 values for each of the examined drugs at 0.5, 1.0, and 2.0 is given in Table 2.2.

### **2.3.2 Electromechanical Window**

Looking at the values, high risk drugs were classifed by 20% decrease compared to control; fecainide, ibutilide, and quinidine returned values below this threshold, indicating high risk, while metronidazole and mexiletine returned values above it. Finally, for ranolazine the decrease was slightly less than 20% which, in combination with the remaining biomarkers, led us to classify it as an intermediate risk drug. For all the drugs, the TdP risk classifcation based on the was the same as 90 risk classifcations. In four out of the six analyzed drugs, the pro-

Fig. 2.3: Action potentials for control and six tested drugs at 1.0. Based on these results, ibutilide, quinidine and fecainide are classifed as high risk. Ranolazine is classifed as intermediate risk for this drug concentration.

gressively decreased with higher drug concentration compared to control. However, addition of metronidazole at 0.5 increased the initially while for higher concentrations of metronidazole (2.0) the decreased again. The addition of mexiletine changed the by only 3 ms while both the 90 and 90 remained similar to control. All results are shown in Table 2.2.

### **2.3.3 Maximum Active Tension**

Applying metronidazole, the maximum active tension decreased from 55.0 kPa to 28.8 kPa, 22.9 kPa, and 17.8 kPa for 0.5, 1.0 and 2.0, respectively. This represents a decrease in , of 48%, 58%, and 68%. For ranolazine, the , increased by 45% at 1.0 compared to control. In the remaining drugs, the maximum active tension remained similar or slightly increased compared to control.

Fig. 2.4: Action potential for control and quinidine at 0.5, 1.0, and 2.0. An increase in PC highlights the drug concentration efect on the action potential duration.

### **2.3.4** *In Silico* **and** *in Vitro* **TdP Risk Classifcations**

We classifed all six drugs using the evaluation process described in section 2.2.3. Flecainide, ibutilide, and quinidine were determined to be high risk both in the previous classifcation and in our electromechanical model (based on 90 and values). Based on the maximum active tension, the risk of metronidazole was changed from intermediate (CiPA indication) to high risk. On the other hand, following the 90 and values it would have been classifed as low risk. Mexiletine was classifed as low in the CiPA as well as in our study. Lastly, ranolazine was previously predicted as low risk but was classifed as intermediate risk using our model based on both the 90 and values. Even though the increase in , for ranolazine was the largest of the analyzed drugs, potentially having an efect on contraction and tissue stress, increases in a,max were not considered in TdP risk classifcation.

### **2.4 Conclusions**

In this study we analyzed six drugs, fecainide, ibutilide, metronidazole, mexiletine, quinidine, and ranolazine. Three of these drugs were classifed as high risk, one as intermediate risk, and two as low risk according to CiPA. Using our fully coupled electromechanics model, four drugs were classifed as high risk, one as intermediate risk and one as low risk. An overview of this distribution is given in Table 2.3. After performing our analysis, we concluded that some drugs may belong to diferent risk categories depending on the individual parameters considered. For example, ranolazine was previously classifed as low risk; in our simulation, both the 90 and biomarker values for 1.0 fell in the range of intermediate risk, while for 2.0 it was predicted high risk. In general, for all six drugs, both the 90 and the biomarkers indicated the same risk categories. However, the risk classifcation based on , difered from the other biomarkers. For metronidazole, a low TdP risk was observed based on the electrophysiological biomarkers. Despite that, there was a decrease of 58% in , which can lead to severe problems in heart function.

Table 2.3: Predicted versus *in vitro* risk classifcations for the analyzed drugs. Out of the six tested drugs, two drugs were classifed diferently based on the electromechanical results.

We conclude that the addition of a mechanical biomarker, the maximum active tension, is valuable in the prediction of cardiac risk stratifcation. Additionally, using the electrophysiological biomarkers, SimCardEMS succeeded in predicting the same TdP risk category as previously determined for four of the six drugs.

2 Electromechanical Testing in TdP 29

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 3** *In silico* **Investigation of Sex-Specifc Osteoarthritis in Human Articular Chondrocytes**

Khoa Ngo, Nathaniel T Herrera, Milda Folkmanaite, Kei Yamamoto, and Mary M Maleckar

**Abstract** Osteoarthritis (OA), a progressive degenerative disease of cartilage in joints, is the most common cause of chronic disability in older adults. While OA is mostly considered an age-related pathology, women have a 1.5-fold higher risk of developing OA relative to men and experience more severe symptoms. Yet, they remain underrepresented in musculoskeletal research and clinical trials. Responsible for cartilage formation, articular chondrocytes experience physiological changes in OA, but the functional implications of such alterations remain largely unexplored due to difculties in acquiring the data experimentally. Through reparameterization, we expand a mathematical chondrocyte model to investigate sex-specifc OA pathogenesis. We performed sensitivity analysis to address the impact of ion channel activity in healthy and OA chondrocyte populations. Simulations show that in healthy female chondrocytes, the resting membrane potential is more depolarized than in healthy male chondrocytes, suggesting potential sex-specifc emergent physiological diferences in articular chondrocytes. In both sexes, the resting membrane potential of healthy chondrocytes is most sensitive to − , −, and −, but in OA it depolarizes and becomes sensitive to , and −. Developed and evaluated against experimental data, our articular chondrocyte OA electrophysiolog-

Khoa Ngo Department of Physiology and Membrane Biology, University of California, Davis, USA

Nathaniel T Herrera Department of Pharmacology, University of California, Davis, USA

Milda Folkmanaite Department of Physiology, Anatomy and Genetics, University of Oxford, UK

Kei Yamamoto Department of Computational Physiology, Simula Research Laboratory, NO

Mary M Maleckar (corresponding author)

Department of Computational Physiology, Simula Research Laboratory, NO e-mail: mmaleck@ simula.no

ical model can be used to further study OA pathology and sex-specifc pathological OA changes.

### **3.1 Introduction**

Mammalian chondrocytes are found within intervertebral discs and articular cartilage, where they maintain the extracellular matrix (ECM) and produce the cartilage matrix [1]. The ECM mainly consists of collagen fbers that transmit force and dissipate energy, proteoglycans that withstand compressional forces, and synovial fuid to reduce friction between bones. In normal function, healthy chondrocytes respond to outside stimuli and damaged articular cartilage, where they increase production of specifc ECM complexes [2]. Functional, mature, healthy articular chondrocytes are necessary for the maintenance of cartilage and therefore, physiological joint motion [3]. Chondrocytes function to maintain a homeostatic environment in the articular cartilage. Pathophysiologic conditions can lead to the development of osteoarthritis (OA); however, the mechanism by which this occurs is poorly understood. Osteoarthritis is the most common form of arthritis and is a disease characterized by the degeneration of cartilage and underlying bone [4]. Women are 1.5x more likely to develop OA, and experience more severe symptoms as compared to men [5, 6, 7]. In OA, chondrocytes secrete increased levels of infammatory cytokines, actively produce proteoglycans and collagen type II to recover the degeneration of the ECM, and become hypertrophic [8]. The present study on human, sex-specifc OA is modeled after the specifc cell physiology of the human articular chondrocyte.

Recent research suggests that maintaining a stable resting membrane potential is essential for chondrocytes in articular cartilage to be able to withstand pressure and force changes [9]. Despite chondrocytes being non-excitable cells, they undergo robust ion-channel mediated changes in response to OA. This is further supported by experiments where ion-channels contributing to the resting membrane potential were blocked, resulting in a decrease in the production of matrix mRNAs, proteins and glycosaminoglycans [10, 11]. Additionally, channel blockers inhibit chondrocyte proliferation and increase apoptosis [12, 13]. During OA, chondrocytes have reduced capability to respond to changes in the extracellular environment, specifcally with respect to ion transport; this is recognized by i.e. a defciency in volume regulation [14]. Chondrocyte diameters vary between 7 and 30 µm, making technical electrophysiological experimental studies difcult [15]. Therefore, there is limited experimental data investigating the electrophysiological diferences of chondrocytes in normal physiology and pathophysiology.

Therefore, to gain a better understanding of the human chondrocyte channelome in the context of OA, we have developed an OA chondrocyte model by further expanding a previously established mathematical model of the human articular chondrocyte [3, 16]. Mathematical models can play a key role in the dearth of physical data and in hypothesis generation, among other utilities. In the absence of constraining data, models can be used to generate valuable hypotheses about the underlying processes at work, which can then be tested in carefully designed experiments and/or against data from other sources. This initial model of the human articular chondrocyte features various ion channels describing the resting membrane potential and ion movement across the cell membrane. Previous iterations of the model focused on potassium fuxes across the cell membrane, and included functions specifcally for the integration of the sodium/potassium pump, an active transport pump that establishes an electrochemical gradient across the cell membrane. In this study, the model is expanded to include sex-specifc OA by implementing the diferences in ion channel conductances observed between males and females, as well as alterations of ion channel expression in OA disease. Populations of models using a log-normal distribution were also implemented to accurately represent biological variability, and subsequent parameter sensitivity analysis was performed to identify the currents that infuence chondrocyte resting membrane potential to the greatest degree in both health and disease.

In summary, we have developed a mathematical OA chondrocyte model to gain a better understanding of OA pathogenesis through exploratory simulations, driving forward hypothesis development in the context of relative scarcity of human chondrocyte electrophysiological data across the channelome. As OA is more prevalent in females than males, we also implemented sex diferences in our updated model. Our main goals in the present study were to: (i) update and scale the baseline male model (1) to generate a (2) disease OA-male model, (3) a control female model, and (4) a disease OA-female model, (ii) utilize a population of models approach to determine which ion channels were the main contributors to resting membrane potential in human articular chondrocytes during normal function and in OA.

### **3.2 Methods**

To study the efects of OA on the resting membrane potential, we employed the chondrocyte model and simulation protocols developed by Maleckar *et al.* [3] and further improved by Fischer-Holzhausen *et al.* [17]. The model components are shown in Figure 3.1. has been scaled to be within 0.1-0.6 of the original value, corresponding roughly to 0.9–5.35 pA/pF, as prior work has shown this range is more likely to present a physiologically-relevant model regime [17].

Following the Hodgkin-Huxley formalism, the cell membrane is modeled as a capacitor coupled in parallel with ion channels represented by resistors. The timeevolution of the transmembrane potential is given by the summation of all electrogenic transport processes.

$$C\_m \frac{dV}{dt} = -(I\_{KDR} + I\_{NaK} + I\_{NaCa} + I\_{Ca-ATP} + I\_{K-ATP})\tag{3.1}$$

$$+ I\_{K2pore} + I\_{Nab} + I\_{Kb} + I\_{Clb} + I\_{KCa})\,,$$

Fig. 3.1: Schematic illustration of the model adapted from Maleckar *et al.* [17], showing the principal ion-selective channels, exchangers, and pumps expressed in human chondrocytes.

where is the chondrocyte capacitance (6.3 pF). A set of three transport equations gives the time-derivatives of the intracellular ionic concentrations of + , + , and 2+

$$\frac{d\left[Na^{+}\right]\_{i}}{dt} = -\frac{I\_{Nab} + 3I\_{NaK} + 3I\_{NaCa} - I\_{NaH} + I\_{Clb}}{vol\_{i} \times F},\tag{3.2}$$

$$\frac{d\left[K^{+}\right]\_{i}}{dt} = -\frac{I\_{Kb} - 2I\_{NaK} + I\_{KDR} + I\_{K2pore} + I\_{KCa} + I\_{K-ATP}}{vol\_{i} \times F},\tag{3.3}$$

$$\frac{d\left[Ca^{2+}\right]\_i}{dt} = -\frac{I\_{Ca-ATP} - 2I\_{NaCl}}{2vol\_i \times F} - 0.045 \times \frac{dO\_c}{dt} \,, \tag{3.4}$$

where is the internal volume of the chondrocyte (0.005884 mL), is the Faraday constant (96,485 C/mol), and is the fraction of intracellular calmodulin bound to 2+ .  For example, the NaK exchangers exchange three + ions for two + ions.

### **3.2.1 Modeling Ionic Changes Induced by Osteoarthritis in Male and Female Chondrocytes**

To gain mechanistic insight into the pathogenesis of OA, we implemented changes in ion channel expressions as measured by Karlsson *et al.* [18]. We further expanded our investigation by integrating sex-specifc diferences in ion channel activity into our model. Due to the lack of experimental data on female chondrocytes, we compiled a list of sex-specifc diferences in ion channel expression based on mRNA and ion channel subunit expression data collected from epicardial cardiomyocytes. In our model, we modeled the efect of changes in ion channel expression on their conductances (), as displayed in Table 7.2.

Table 3.1: Sex-specifc subcellular ion channel conductances relative to the male baseline model

Implemented as change to the Q10 parameter which defnes the efect of temperature on ion channel kinetics [19, 20].

### **3.2.2 Generation of a Population of Models**

To capture biological variability and investigate parameter sensitivity, we used the published chondrocyte model [3, 17] as a baseline to build our population of 1,000 models by randomly modifying parameters corresponding to conductances of ion current and maximal rates of ion transports. A relatively large sample size of 1,000, typically used in similar population-based studies [21, 22], was chosen to allow even small efects, often indistinguishable in small sample sizes, to reach statistical signifcance. Prior studies have shown that biological parameters, when randomly selected, are distributed according to log-normal frequency curves instead of normal ones [23]. For each channel in the model, we varied its conductance by sampling from a log-normal distribution centered around their respective baseline value with the standard deviation set to 0.15.

### **3.2.3 Model Parameter Sensitivity Analysis**

We performed parameter sensitivity analysis on the generated model population to assess the sensitivity of the model output (e.g., the resting membrane potential) in response to channel conductance variations. The degree to which perturbing a set parameter can infuence the results are quantifed as a set of regression coefcients, which are obtained by performing multivariable partial least squares regression on the dataset. Partial least squares regression was chosen over standard multivariable regression as many independent parameters (e.g., ion channel conductances) were used to predict a smaller set of dependent variables (e.g., the resting membrane voltage) [24]. This methodology was frst used in cardiac electrophysiology by Sobie *et al.* [25] and has been widely used in other felds of biology [21, 26].

### **3.3 Results**

### **3.3.1 Modeling the Impact of Sex-Specifc OA on Resting Membrane Potential**

In Figure 3.2, we implemented the experimentally-observed changes in ion channel expression in control and OA chondrocytes, accounting for sex diferences as described in Table 7.2. Figure 3.2 shows that the healthy female chondrocyte maintains a more depolarized resting membrane potential (-58.21 ± 5.75 mV) as compared to the healthy male chondrocyte (-69.11 ± 4.71 mV). Electrophysiological diferences instigated by OA result in depolarization in both male (from -69.11 ± 4.71 mV to -53.87 ± 3.21 mV) and female (from -58.21 ± 5.75 mV to -49.03 ± 1.34 mV) chondrocytes. This overall depolarizing efect of OA remodeling on the resting membrane potential has previously been reported in experiments [27].

Fig. 3.2: Resting membrane potential in control and OA chondrocytes for (A) male and (B) female models. Changes in ion channel expression in OA induce resting membrane depolarization.

### **3.3.2 Population of Models**

Shown in Figure 3.3, and as outlined previously, we further expanded our analysis to investigate population variability by generating 1,000 models for each condition with channel conductances varied according to a log-normal distribution. For all 1,000 parameter sets, the steady-state solutions were reached within the simulation period of 50,000 seconds. Figure 3.3 shows that in male and female OA chondrocytes, intracellular sodium quickly becomes depleted and remains very small, close to zero, for the rest of the simulation. We calculated the mean resting membrane potentials for all conditions and compared them against the reported experimental values as measured in experiments, shown in Figure 3.4. In all cases, the membrane potentials of OA chondrocytes are relatively more depolarized compared to their healthy counterparts. Comparing male and female chondrocytes, although the control resting membrane potential values slightly difer (females are more depolarized compared to males), OA male and female chondrocytes depolarize to a similar value (approximately -51 mV).

### **3.3.3 Parameter Sensitivity Analysis**

Sensitivity analyses were performed on the chondrocyte populations to reveal how ionic modulations induced by OA might afect the sensitivities of underlying electrophysiological properties in the models. Shown in Figure 3.5, in the healthy male chondrocyte population, , − , −, and − are predicted to have the most critical impact on the resting membrane potential. In this case, − is the strongest depolarizing current and − is the strongest hyperpolarizing current. In male chondrocytes expressing OA, only and − remain the most infuential currents. Compared to the healthy case, the impact of − and − are severely reduced in OA, while gains a 1.6-fold depolarizing infuence.

The healthy female chondrocyte population shares some similarities with their male counterparts, as the resting membrane potential is still the most sensitive to , − , −, and −. However, − now plays a slightly more dominant role as a hyperpolarizing current. Compared to male chondrocytes, gains a more signifcant role in establishing the resting membrane potential than in males (2.5-fold increase in importance). is the strongest depolarizing current and − is the strongest hyperpolarizing current. Surprisingly, in female chondrocytes expressing OA, the majority of the currents no longer have as much of an impact on the resting membrane potential. Nonetheless, and − each still retain a minor infuence.

Fig. 3.3: Membrane voltage and internal ion concentrations for (A) male and (B) female population.

### **3.3.4 Inhibiting Restores Normal Resting Membrane Potential in OA Chondrocytes**

Per sensitivity analysis results, has the largest impact on depolarizing the resting membrane potential (i.e., having positive regression coefcients) in both male and female chondrocytes in OA. To restore the membrane potential back to near control values, we applied current block treatment then simulated the membrane potential in male and female OA populations of 1,000 models each. Figure 3.6 shows that applying 45% block to male and 55% block to female OA chondrocytes restored the mean resting membrane potential of the population back to the control level. Notably, previous clinical studies have revealed that

Fig. 3.4: Mean resting membrane potential in control and OA chondrocytes for male and female populations (1,000 models each) compared against experimental data [27]. To account for diferences in data collection methodology, experimental resting membrane potential values were scaled to match the mean of the male population.

Fig. 3.5: Parameter sensitivity analysis on (A) male and (B) female population of sex-specifc chondrocyte models. The regression coefcients measure the relative impact of a particular ionic current on the resting membrane potential. A positive regression coefcient indicates the associated current is a depolarizing current, while a negative regression coefcient indicates a hyperpolarizing current.

blockers Ouabain and Digoxin protected against OA and relieved OA-associated pain [28].

Fig. 3.6: Resting membrane potential in control, OA, and OA chondrocytes with partial block of for male and female model population (1,000 models each). Blocking (45% for male and 55% for female population) helped restore the resting membrane potential back to normal levels in OA chondrocytes.

### **3.4 Discussion**

OA, a debilitating disorder afecting 1 in 7 individuals in the United States in their lifetime [29], has no current curative treatment. The prevalence, as well as disease severity, is more pronounced in female patients [30]. It is known that chondrocytes play an essential role in maintaining healthy cartilage, and that their function is impaired in disease [31]. As experimental investigation of human chondrocytes is challenging due to the lack of human control samples, we have further expanded a previously-developed mathematical model of a chondrocyte channelome by introducing experimental changes observed in ion channel electrophysiology, as well as temperature and capacitance [3]. The updated model is user-friendly and freely accessible (github.com/k-ngo/Chondrocyte).

Simulations revealed that the chondrocyte membrane becomes depolarized in OA, which agrees with experimental data [32]. This suggests that the changes in ion channel expression together with temperature and capacitance are sufcient to capture the main pathological alterations observed in OA chondrocytes related to their membrane potential. In the future, it remains important to further validate the model experimentally, as well as to investigate the primary cause of the alterations in ion channel expression observed in OA chondrocytes. These are likely to involve infammatory processes associated with OA [33].

Comparisons between male and female chondrocyte models reveal intriguing phenotypic diferences. First, the resting membrane potential is more depolarized in females compared to males in controls. The main diferences between male and female healthy chondrocyte models include reduced (by 25%), − (by 95.2%), (by 2.8%) and (by 5.7%). The investigation of the main factors important for maintaining membrane potential visualized in Figure 3.5 revealed that, in the present model, calcium ATPase and background sodium and potassium currents contribute the most for maintaining healthy resting membrane potential in males, with an addition of in females. It could be that the relative depolarization of the female chondrocyte membrane driven by the aforementioned changes in the model predisposes female patients to OA, although additional experimental data is needed to further investigate this hypothesis.

OA chondrocyte model simulations show that resting membrane potential is depolarized in disease, following experimental data [32]. This efect is seen in both male and female OA models. Our results also suggest that chondrocyte resting membrane potential in female patients could be more depolarized compared to that in males, which could lead to higher depolarization in OA (mean resting membrane potential is -49.03 ± 1.34 mV in the female and -53.87 ± 3.21 mV in the male OA model). The resting membrane potential in OA was around 17.5% more depolarized for female and 24.5% for male model compared to respective healthy resting membrane potentials. In our OA model, current plays a prominent role in regulating membrane voltage: in both males and females with OA, increased becomes the most important contributing factor to the overall depolarization of the resting membrane potential, as revealed by sensitivity analysis. Inhibiting might therefore provide a useful approach for reversing resting membrane alterations in OA and, more broadly, could be a useful target in modulating disease pathogenesis in OA.

 inhibitors are well characterized for e.g. cardiac disease treatment, and new selective inhibitor drugs are currently in the pipeline [28]. To test whether block could be useful in overcoming OA-associated changes in resting membrane potential, we simulated block. With varying degrees of inhibition, the resting membrane potential in disease can be restored to healthy, pre-OA levels in both males and females. Therefore, re-purposing blockers feasibly provide a useful new therapeutic approach for OA treatment. blockers have, in fact, already been applied for OA treatment: clinical studies have reported positive outcomes in patients, revealing that blockers, Ouabain and Digoxin, induced a reduction in pain and elicited an apparent OA-protective efect [28], supporting our simulation results implying blockers as a potential treatment for OA.

Given this putative critical role of Na+ concentration in healthy and OA chondrocytes, novel future experiments should consider thoroughly probing Na+ balance in human chondrocytes. For instance, Figure 3.3 shows that in male and female OA chondrocytes, intracellular sodium quickly becomes depleted and remains very small, close to zero, for the rest of the simulation, suggesting that alternate sodium regulation not currently accounted for in base models may apply in OA. Model development should additionally move towards incorporating new channel species to regulate Na+ concentrations in control and in OA, including epithelial sodium channels (ENaC), which are not currently included in the model [31, 34]. Other work has also focused on voltage-gated sodium channels in chondrocytes, which may also be considered.

Naturally, the lack of available biological data from new electrophysiological and other experimental studies presents a challenge for the validation of both male and female models, both in control and in OA. Further experimental studies will be instrumental in determining the ranges of parameters for the model, as well as for simulated chondrocyte behavior evaluation. In the future, further experimental data will be helpful in further developing and expanding the current chondrocyte model. While the current study, synthetic in nature, limits its direct translation to clinical application, the present work retains utility in terms of physiological exploration and hypothesis generation.

In summary, further incorporating additional species key in the chondrocyte channelome into the model and incorporating new data from emerging experiments will permit even greater insights into the broad and nuanced role of membrane dynamics in chondrocyte function and pathophysiology.

### **3.5 Conclusions**

We have further expanded a previously published chondrocyte mathematical model enabling simulations of the osteoarthritic chondrocyte channelome. Implemented changes were based on published data, and model simulation results show that the resting membrane potential in chondrocytes becomes depolarized in OA in both male and female chondrocytes, which agrees with the experimental data acquired from males. Our sex-specifc chondrocyte model reveals diferences in resting membrane potential between males and females that could potentially contribute to the higher prevalence of OA in female patients. Finally, sensitivity analyses revealed the main currents responsible for maintaining resting membrane potential in chondrocytes and potential novel therapeutic targets for OA treatment. Our OA chondrocyte electrophysiological model provides an accessible tool for subsequent studies of OA pathogenesis and drug targeting.

3 Sex-specifc OA in Chondrocytes 43

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 4 Recapitulating Functional Heterogeneity in Electrophysiologically Active Tissues**

Meye Bloothooft, Joseph G Shuttleworth, Gabriel Neiman, Ishan Goswami, and Andrew G Edwards

**Abstract** Inter-cellular heterogeneity is central to the dynamic range and robustness of function in many tissues, particularly electrically excitable tissues. In pancreatic islet -cells, inter-cellular heterogeneity underlies the range of insulin response to glucose. In human-induced-pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs), inter-cellular heterogeneity presents a key challenge for drug screening applications. In this study, we assess the ability to reconstruct inter-cellular heterogeneity *in silico* by applying a "population of models" (PoMs) framework, i.e. collections of computational cells created via Monte Carlo variation of model parameters. We defne parameter variation based on experimentally observed heterogeneity in properties such as ion current conductances and enzymatic afnities. We then assess the accuracy of those reconstructions, based on the degree to which variation in PoM outputs (e.g. action potential duration) matches experimentally observed variation. We report that this "ground-up" approach underestimates functional heterogeneity in the hiPSC-CM population, but overestimates it in adult human cardiomyocytes. In contrast, the -cell PoM captures three distinct and physiologically relevant subclasses of -cell function. In the future, we expect PoM approaches like these will

Meye Bloothooft

Department of Medical Physiology, Division of Heart and Lungs, University Medical Centre Utrecht, NL

Joseph G Shuttleworth Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, Nottingham, UK

Gabriel Neiman Department of Bioengineering, University of California, Berkeley, USA

Ishan Goswami

Department of Bioengineering, University of California, Berkeley, USA, Department of Materials Science and Engineering, University of California, Berkeley, USA

Andrew G Edwards (corresponding author)

Department of Computational Physiology, Simula Research Laboratory, NO e-mail: andy@ simula.no

permit incorporation of realistic cellular heterogeneity in detailed models of intact tissues, and thereby aid development of sophisticated tissue-engineered platforms for therapeutics.

### **4.1 Introduction**

Functional heterogeneity among cells within a tissue is a characteristic and relevant feature of living systems. Even within a genetically homogeneous population, individual cells often show a considerable amount of phenotypic heterogeneity [1], and this feature has implications in the plasticity and robustness of the organisms. Plasticity allows adaptive cellular responses to perturbations in the environment while maintaining a robustness in functional behaviour of the tissue even after a perturbation or insult. For example, when a region of the heart becomes infarcted (dies), non-infarcted healthy areas will adapt and increase their force of contraction to maintain blood fow to the body. This ability to adapt is a form of robustness of organ function, but must be achieved in a manner that maintains the heart's electrophysiologic stability. Similarly, in obesity induced diabetes, hyperlipidemia and hyperglycemia exert excessive insulin demand on pancreatic islet -cells. In some -cells this prolonged exertion leads to exhaustion of insulin reserve and cell death. Thus, increased demand for insulin must be achieved via an adaptive response in the remaining viable -cells, which are phenotypically distinct. In both of these examples, cellular heterogeneity is central to the plasticity and robustness that allows both organs to maintain function after pathophysiological challenge. Functional heterogeneity may also broaden the range of responses to drugs [2]. For example, glibenclamide has diferential efects on pancreatic islets owing to functional heterogeneity of their -cell populations [3]. Thus, the study of cellular heterogeneity for both basic understanding of living systems and discovery of therapeutics has motivated the development of experimental techniques and the construction of computational models to simulate the dynamics and efects of cellular heterogeneity [4].

A range of single-cell experimental techniques are available to investigate this heterogeneity. These range from transcriptome assessment via single cell RNA sequencing, to classical cell electrophysiology via patch clamp. While these approaches are generally both time and labor intensive, data collected by those methods can be stored and interrogated via efcient modelling frameworks, such as the population of models approach (PoM) [5, 6]. These approaches simulate the dynamics and functional implications of heterogeneity in cell populations, including responses to pathologic insults and drug challenges.

In the PoM approach, heterogeneity is implemented by varying parameters (such as the ionic conductance of a single or multiple ion channel/s) of a baseline model, which itself captures the average characteristics of a particular cell type. Each unique confguration of parameters creates a new instance of the cell model, which is analogous to a new "member" cell of the population. In general, this procedure has been performed by frst specifying equal distributions over which all parameters of interest (e.g. all ionic conductances) can vary, and then selecting only members of that population that exhibit behaviour within the experimentally observable range [5, 6]. For example, a cardiac myocyte population could be constructed by frst allowing all ionic conductance parameters to be sampled from Gaussian distributions with an arbitrary standard deviation of 30% of each mean parameter value in the baseline model - this sampling strategy involves homogeneous variability across parameters. Some members of this PoM will exhibit parameter combinations that are biologically unrealistic, and thus possibly exhibit behaviour outside the experimentally observable range. To eliminate those models, it is common to reject models that produce outputs or features, such as a cardiomyocyte action potential duration, that are outside the range of experimental observations for the respective features. Continuing the example, a cardiomyocyte population member will be deemed non-viable if the action potential duration is far longer or shorter than observed experimentally. Similarly, a non-fring -cell model will be rejected from the islet PoM. Because, in these examples, the PoMs are generated from arbitrary and homogeneous parameter distributions, they will exhibit appropriate behavioural characteristics, but may not involve truly biologically representative parameter combinations. Thus, the behaviours may not result from truly representative underlying dynamics.

In this work, we apply a diferent approach, with the intention of introducing more biologically realistic variability in PoMs of three diferent cell types: human adult cardiomyocytes, hiPSC-CMs and primary islet -cells. By introducing variability to the input parameters that is grounded in *experimentally observed variability for each parameter*, we hope to span the biologically reachable confgurations and limit the number of unrealistic confgurations that are generated. That is, we use published experimental variation (standard deviations) for ionic conductance parameters and enzymatic activity/afnity parameters in order to apply parameter-specifc heterogeneity rather than variation that is homogeneous across the varied parameters. Based on these constraining conditions, we demonstrate the feasibility of PoMs to generate a heterogeneous population representing diferent functional phenotypes.

### **4.2 Methods**

An illustrative example of our approach in generating PoMs is shown in Figure 4.1. To recapitulate the functional heterogeneity in the tissue PoMs, we introduced parameter-specifc variability into our models. To do this, we identifed a set of parameters which we expect to account for much of the inherent biological variability infuencing electrophysiologic function in these cell types. These parameters are a selection of maximal conductance (or fux) parameters for both the cardiomyocyte models (hiPSC-CM and adult) and the -cell model. Additionally for the -cell, a selection of kinetic parameters for enzyme reactions in the glycolysis cycle were varied to introduce metabolic heterogeneity. We chose to endow each of these pa-

Fig. 4.1: Illustrative example of generating population of model (PoM) of human cardiomyocytes and islet -cells. Experimental data obtained from literature was used to generate log-normal distributions of electrophysiological and/or metabolic parameters used as inputs to computational model. Sets of features are obtained from the simulations that are then used to generate viable candidates for the PoMs representative of cardiomyocyte or -cells via constraining conditions.

rameters, , with a log-normal distribution, such that log ∼ ¯ (0,2 ) where ¯ is the original, unmodifed parameter from the baseline model and is a standard deviation derived from the literature. We term the random variable ¯ the *scaling factor*. To assign values for the various , we surveyed the literature for measured variability of each parameter of interest. To mitigate the impact of systematic diferences among the means of diferent datasets, and to keep our parameter distributions centred on the original parameter values, we computed the coefcient of variation for each parameter ():

$$
\sigma\_p \coloneqq \frac{\sigma\_{\text{lit}}}{\mu\_{\text{lit}}} \tag{4.1}
$$

where lit and lit are the mean and standard deviation, respectively, of the observations found in the literature.

The published mathematical models chosen as the baselines for the adult-CM and the hiPSC-CM were developed by [7] and [8], respectively. For the -cell PoM, the electrophysiological model developed in [9] and glycolysis model developed in [10] were used to generate the heterogeneous population of -cells. In the following subsections, we will detail the generation of parametric distributions and the constraining conditions via feature extraction for each PoM.

### **4.2.1 Cardiomyocytes**

The O'Hara-Rudy model (ORd, [7]) is a widely-used mathematical description of the average human left ventricular cardiomyocyte. We specifcally used the epicardial version for our adult human baseline model. Likewise, the baseline model for hiPSC-CMs was that published in [8]. Both models use the following fundamental equation for calculation of time-varying membrane potential:

$$C\_m \frac{dV}{dt} = -\sum\_i I\_i \tag{4.2}$$

where is the capacitance of the cell membrane and each is a current describing the transmembrane current carried by each species of ion channel present in the model. For each ion channel, there is a maximal conductance parameter , or analogous permeability parameter, such that = ˆ for some ˆ . ˆ is in turn described by a system of ordinary diferential equations that specify the voltage, ion concentration and time-dependencies of ion channel function. For all currents, we only vary the maximal conductance (), and not the parameters of the remaining equations contributing to ˆ . The formulation for sarcoplasmic reticulum (SR) calcium reuptake is diferent to those of the ion channels, but an analogous term that refects the maximal rate of SR calcium reuptake ¯ SERCA. As for the ion channel conductances, this parameter was varied to introduce heterogeneity in this key element of intracellular calcium handling.

Table 4.1: **Scaling parameters for the cardiomyocyte PoM**

Sources to construct coefcient of variation () listed by the reference numbers

[11], 2 [12], 3 [13], 4 [14], 5 [15], 6 [16] [17], 8 [18], 9 [19], 10[20], 11[21], 12[22], 13[23], 14[24], 15[25], 16[26], 17[27], 18[28], 19[29], 20[30], 21[31], 22[32]

Fig. 4.2: A. The resultant log-normal distributions for the Ks scaling factor, along with histograms resulting from 100 samples, are shown for the representative parametric variation in adult CM and the hiPSC-CM inputs for the models. B. The scaling factors that were sampled in order to create the PoM are shown for two cardiomyocyte models. Parameters from the models that were excluded from the POM are shown in grey, whereas the parameters which were calibrated into the PoMs are highlighted in orange.

Table 4.1 provides the values for each varied parameter for both the adult human myocyte model and the hiPSC-CM model. The variability presented in this table may not completely capture the variability in each maximal conductance parameter. Nevertheless, these values provide a reasonable approximation, and yield plausible distributions for each of the relevant parameters. We took 500 samples from each of these distributions and the results are shown in Figure 4.2A alongside the relevant sampling distributions.

The PoMs generated by varying the input parameters as summarized in Table 4.1 were solved via MATLAB stif-solver ode15s, and allowed to equilibrate by simulating for 500s before computing the output features. Calibrated PoMs for both the adult CM and hiPSC-CM were obtained by rejecting models that produced features outside the experimentally observed range. We identifed a collection of features from the literature: *maximum upstroke velocity (voltage)*, *minimum diastolic potential*, *voltage amplitude*, *action potential duration (APD*90*)*, *maximal departure velocity for Ca2+*, *diastolic Ca2+ concentration*, *Ca2+ transient amplitude*, *Ca2+ transient time-to-peak*, and *time constant of Ca2+ decay*. Representative scaling factors sampled from the output of the models obtained via the constraining criterion are shown in Figure 4.2B. In the results section, we discuss the typical features obtained and constraining criteria to generate calibrated PoMs for the cardiomyocytes.

### **4.2.2 Islet -cells**

A robust secretion of insulin by the -cell upon glucose challenge relies heavily on the coupling of the metabolic oscillations in the glycolysis pathway with electrical oscillations (e.g. isolated action potentials and action potential bursting). Recognizing this, the -cell PoM was generated via a coupling of the electrophysiological model developed in [9] and glycolysis model developed in [10]. A coupling of these models was demonstrated in [33], albeit using single mean values for the electrical and glycolytic component parameters. A simplistic illustration of this coupled model is shown in Figure 4.3. To delineate contributions made by the glycolytic and electrophysiologic systems to the overall functional heterogeneity in our -cell populations, we divide the parametric distributions into two components: ionic and glycolytic.

Fig. 4.3: Schematic representation of a -cell's metabolic and electrophysiological circuit involved during the glucose-stimulated-insulin-secretion. Included in the schematic are representative formulae of the glycolytic component's enzyme rate equations and a few parameters (in red) that are changed to generate PoMs. Lognormal distributions for two of the 14 parameters are shown in the bottom right hand side of the fgure.

Details of the ionic models can be found in the cited articles, but briefy electrical oscillations in the system were modeled as described by Equation 4.2. The maximal conductance () of 8 ion channels were varied based on experimental data in the literature (references found in [34], [35], and [36]). These variations, reported via coefcients of variation () and their sources, are summarized in Table 4.2.

The glycolysis model comprised enzymatic reactions from glucokinase to glyceraldehyde-3-phosphate. The rate of ATP generation through these steps and via downstream oxidative phosphorylation in the mitochondria is introduced in a phenomenological manner via the variable, *a*, which represents ATP level.

$$\frac{da}{dt} = V\_{GAPDH} - k\_A a \tag{4.3}$$

where is the metabolic fux due to GAPDH and is a phenomenological time constant.The ion channel conductance of channel is then altered via:

$$\text{g}\_{KATP} = \frac{\text{\hat{g}}\_{KATP}}{1 + a} \tag{4.4}$$

In our study, we varied the kinetics of the enzymatic reactions by varying the limiting rates , of glucokinase, phosphofructokinase, and glyceraldehyde 3-P dehydrogenase. In addition we varied other kinetic parameters such as the halfactivation concentrations 0.5 . These parameters are highlighted in red in our Figure 4.3. Variations in these glycolytic components of our model and their sources are summarized in Table 4.2. This variability was again assumed to be log-normal distributed and implemented as such for construction of our PoMs. Sources for these datasets can be found in [37], [38],[39],[40].

Since the behaviour of individual -cells is much more diverse than those of cardiac myocytes, the constraining conditions are also less well-established. Thus, we discuss the features used for constraining the -cell PoMs in the results, and we report sub-classes of -cells obtained via our approach. Furthermore, we will analyze the efects of perturbations to the glycolysis and ionic components in maintaining functional heterogeneity.

### **4.3 Results**

### **4.3.1 Cardiomyocyte PoM**

As described in the methods, we allowed each PoM member model to equilibrate for 500s before computing the output features. The last 10 action potentials of each simulation for these models were used to extract features that were compared against the constraining criteria to generate the calibrated PoMs. Constraining criteria were:

#### 4 Functional Heterogeneity 53

### Table 4.2: **Scaling parameters for islet -cell PoM**

Sources to construct coefcient of variation () listed by the reference numbers

[34], 2 [35], 3 [36], 4 [37], 5 [38], 6 [39], 7 [40]

The full populations, showing both the accepted (orange) and discarded (grey) model confgurations are shown in Figure 4.4A, together with the percent of confgurations discarded based on each pairwise combination of constraints in Figure 4.4B. This fgure provides a frst indication of how well variability in the model behaviour is captured by applying a data-defned variation in the model inputs. A large number of confgurations had to be discarded for both models, although the reasons for exclusion were diferent. Specifcally, we retained only 55 of 500 models for our adult ventricular cardiomyocyte PoM, and 60 of 100 models for our hiPSC-CM PoM. For the adult ventricular PoM the vast majority of exclusions were made because the confguration failed either 1 (43% of confgurations) or 9 (40% of confgurations). This can also be seen in Figure 4.4A, where many confgurations fail to repolarize below -65 mV (9), whereas others reach supraphysiologic membrane potentials (either > 70mV or < −100mV). For the hiPSC-CM model the reasons for exclusion were physiologically reachable but not observed in published hiPSC-CM phenotypes. That is, the vast majority of excluded confgurations simultaneously failed criteria 2 - 9, each of which was chosen to address a diferent source of biologically reachable instability. This can also be seen in Figure 4.4A, where most of the discarded (grey models) have attracted to a second equilibrium in membrane potential between −10mV and −20mV. This property of electrophysiologic bistability is a real biological property of cardiac myocytes [41]. It is thus a fair question to ask whether such confgurations should be discarded for the hiPSC-CM PoM, as they may in fact be realistic members of most hiPSC-CM cell populations. Regardless, this initial characterization suggests that the experimentally-defned variation in hiPSC-CM input parameters results in sharp transition to gross electrophysiologic instabilities, most of which remain physiologically plausible. On the other hand, the majority of adult ventricular PoM exclusions resulted from either biologically unreachable phenotypes (V > 70mV or < −100mV), or to an unstable resting potential (minimum V > −65mV).

Fig. 4.4: **A** The action potential and calcium transient phenotypes from the adult CM and hiPSC-CM PoMs. Traces of cells accepted into the calibrated PoMs are in orange, while those discarded are shown in grey. **B** Tables summarizing the % of models excluded by the constraining criteria. Diagonal values show the percentage of models excluded by each criterion in isolation, whereas the of-diagonal entries show the % of confgurations for which each combination of two excluding criteria were present. The meaning of each criterion, i.e. , is explained in the main text.

Figure 4.5 shows the relationships between major features of the fnal calibrated PoMs for the two cardiac models, and provides a basis for comparing systematic diferences in the phenotypic span of each PoM. The stronger relationships for the hiPSC-CM PoM between Ca2+ handling features and AP amplitude indicate the larger role of sarcolemmal fuxes in determining Ca2+ features. In contrast, the adult

Fig. 4.5: Scatter plots of feature pairs from both PoMs. The ftted distributions are shown along the diagonals using a kernel density estimator. Noteworthy diferences between the two PoMs are boxed in the red (adult) and purple (hiPSC-CM).

ventricular PoM Ca2+ dynamics are more internally determined, as demonstrated by stronger relationships boxed in green. These diferences refect intrinsic characteristics of the baseline hiPSC-CM and adult ventricular cardiac myocyte models. They are a well-known diference in the biological underpinnings of the two cell types and refect the generally less "mature" Ca2+ handling phenotype of hiPSC-CMs relative to adult ventricular myocytes.

Finally for the cardiac PoMs, Figure 4.6 shows important diferences between the two calibrated PoMs in terms of degree to which the variability in their AP and Ca2+ handling features refect those reported in the literature. Specifcally, it is clear that the adult ventricular PoM features (top row, orange histograms) are generally more variable than is reported experimentally (dashed lines). Only the APD90 feature for this PoM appears to approximate reported experimental variation. In contrast, the hiPSC-CM feature distributions (bottom row, orange histograms) are less variable than published reports of these phenotypes (dashed lines). This presents a fundamental distinction, and suggests that literature reports of variation in the underlying currents and fuxes overestimate the true physiologic variability in human adult cardiac ventricular myocytes. In contrast, the same techniques applied to hiPSC-CMs underestimate observable phenotypic variability. We explore this important observation in the discussion.

CaT amplitude / umolL 1 APD90% / s Fig. 4.6: Histograms showing the variability in the features of the PoM. The frst row shows the adult ventricular cardiomyocyte PoM, and the second row shows the hiPSC-CM PoM. In each plot, the probability density function of a representative Gaussian distribution shows the variation typically found in measurements of these features in the literature. These distributions have been normalised such that they are centered on the means of each histogram.

### **4.3.2 Islet -cell PoM**

The individual models for the -cell PoM were simulated to capture cell behaviour over 10 minutes - the acute phase of glucose-stimulated insulin release. The parameter variations provided in Table 4.2 were used to generate 1000 cells to assess the importance of 3 sources of heterogeneity, by varying model parameters associated with: (a) only glycolytic (enzyme) model components, (b) only ionic (ion channel and transporter) components, (c) both ionic and glycolytic components. In all three cases, each model confguration was simulated under high glucose (10 mM) and low glucose (2 mM) conditions. Thus, the total number of simulations across all 3 cases (3000 cell model confgurations), at both glucose concentrations, was 6000.

Figure 4.7A shows V and cytosolic calcium for representative simulations at high glucose. We see cells 1 through 3 demonstrate similar features as experimentally observable sub-classes of -cells, i.e. spiking, bursting, and plateau phenotypes. Nonfring cells such as those shown in cells 4 and 5 were also observed. As for the cardiac PoMs, in order to generate a calibrated -cell PoM, we have to invoke and apply a set of constraining criteria to reject confgurations that produce features outside the ranges of experimental observations. However, because the -cell phenotype is the most heterogeneous of all cell types assessed here, there exists far fewer such objective constraints in the -cell literature. However, some benchmarks are available. For example, [42] suggested two metrics of activity based on -cell V

Fig. 4.7: Representative simulation outputs and categorizations. **A.** Representative simulation outputs of membrane potential and cytosolic calcium representing different sub-classes of -cells. 1. Spiking cells 2. Bursting cells 3. Plateau cells 4 & 5. Non-fring cells. Scales on the plots represents time of 240s. **B.** Plot of activity fraction vs. Δpeak for the 5 representative cells. **C.** Number of peaks in the Fast Fourier Transform (FFT) of the membrane potential traces for the 5 representative cells.

(the activity fraction and Δpeak), as means of classifying phenotypic sub-classes. Activity fraction was defned as the fraction of time V was above an arbitrary threshold level. The defnition of Δpeak was set as the diference of the two major components of V-time probability distribution (i.e. the diference between "resting" V and "active" V). As for [42], we did not observe clear delineation of the three sub-classes, nor could we readily diferentiate fring vs. non-fring -cells via these metrics (Figure 4.7B). However, we surmised that Fast Fourier Transform (FFT) of the V signal may allow more clear distinction of these classes. Figure 4.7C shows the number of FFT peaks corresponded with the number of dominant harmonics in the V signal for each representative recording in Figure 4.7A.

Based on this analysis, we set the constraining criteria as follows. Any cell with less than 100 peaks in the V peaks, was considered to be quiescent and thus discarded. For cells that had more than 100 peaks in their voltage signal, FFT was performed on the voltage signal to classify them into spiking, bursting, or plateau cells. Based on our initial analyses, cells with less than 20,000 peaks in the voltage FFT were classifed as plateau cells. Bursting cells were classifed as those having anywhere above 20,000 FFT peaks but less than 40,000 peaks. Finally, cells with FFT peaks greater than 40,000 were classifed as spiking cells.

Fig. 4.8: **A.** Percentage of viable cells based on the chosen constraining conditions for six diferent cases (1000 runs each). **B.** Distribution of the sub-classes of -cell from the viable candidates for each of the six cases. HG: High glucose (10 mM). LG: Low glucose (2 mM). Analyses of the cytosolic calcium traces of each of the sub-classes of -cells in the viable population obtained with high glucose: **C.** mean values, **D.** peak values, **E.** Number of FFT peaks. Statistical diferences between the groups were determined by one-way ANOVA followed by Tukey's HSD post-test.

Based on these constraining criteria, we quantifed the viable cells that were not rejected for both glucose concentrations and all 3 cases of heterogeneity (glycolytic, ionic, both), as shown in Figure 4.8A. The greatest reduction in viable cells was observed for glycolytic heterogeneity only, thus highlighting the narrow metabolic parameter range functional -cells operate under and its importance for achieving robust insulin secretion. Another interesting aspect of these data was the number of active cells at low ("resting") glucose (2 mM). This feature of basal activity is unique to human -cells, and underlies basal insulin release at the lower blood glucose concentrations present in humans. Thus, our PoM approach was able to capture this aspect of human-specifc -cell physiology.

Further analysis of the viable candidates allowed us to assess the fraction of each of the 3 PoMs belonging to the sub-classes of -cell phenotypes at each glucose concentration (Figure 4.8B). The lowest fraction of bursting cells was observed in the low-glucose (LG) cases when both components or only ionic components were varied. Upon altering glycolytic components alone, plateau cells comprised the only viable sub-class, thus suggesting the importance of coupling the ionic and glycolytic elements to generate the full range of functional phenotypes observable in human -cell populations. Mean cytosolic calcium of the viable cells within our PoM were not statistically diferent between spiking and bursting cells, suggestive of less pronounced basal shift of calcium signalling (Figure 4.8C). However, plateau cells exhibited a mean cytosolic calcium signature distinct from the other two classes, while peak cytosolic calcium grouped plateau and bursting phenotypes. (Figure 4.8D). Applying FFT analysis to the calcium signals also showed that the number of harmonics in the signal could distinguish the sub-classes of -cells, thus suggesting that this may be a viable means of distinguishing functional heterogeneity in experiments, for which calcium recordings are considerably easier.

### **4.4 Discussion**

In this manuscript, we report the construction of PoMs for three electrophysiologically active tissues that were calibrated against experimental data from the literature. PoMs have proven useful in investigation of cardiac electrophysiological variability and recent studies have furthered the methodology by explicitly incorporating experimental data into the construction of populations of models [43, 6]. In our CM PoMs, we observed that the fuctuations in the input parameters for the models allowed us to generate cells with features comparable to those observed in both primary and hiPSC-derived CMs, albeit with difering coverage of the range of function observable in experiments. In both cardiac PoMs we assumed the main source of variability is difering expression of ion transporters, and that this variation is responsible for the observable variation in their V and calcium cycling features. In hiPSC-CMs particularly, this variation is likely to be pronounced given their immature ("fetal-like") phenotype and the degree to which this can be matured. Importantly, we observed that experimentally-defned variation in the input parameters of the cardiac models resulted in functional output that difered in its range for the two models. That is, for the adult CM PoM, the range of AP and calcium handling features was broader than observed experimentally, whereas it was narrower than experimental for the hiPSC-CM PoM. This is a fundamental observation about the agreement between two types of data recorded in these cell types. For the adult CM it could be that biological covariance (which we have not applied in PoM construction) among the ionic currents reduces the range of reachable AP and calcium features in real cells, or that the relative scarcity of adult human CM data means that these two classes of data are not yet fully internally consistent. In contrast, the greater experimental variability for hiPSC-CM phenotype features suggests either that experimental ionic current measurements are subjected to overly constraining inclusion criteria. Alternatively, our hiPSC-CM model formalism may somehow be constrained to be overly stable with respect to the phenotypic variation resulting from variation in the input currents. Regardless, these questions are fundamental to understanding how perturbations to either adult or hiPSC CMs should be expected to impact those tissues *in vivo*.

The islet -cell PoMs were generated via alterations in both the parameters involved in glycolysis model components as well as the ionic conductance components. We observed that the number of harmonics in V of the simulated computational cells was enough to distinguish diferent sub-classes of -cells observed experimentally, although we believe the proposed constraining criteria requires further validation via prospective screening of experimentally obtained data sets. Furthermore, we observed that the functional heterogeneity of the -cells is dependent on the coupling of the metabolic and ionic components, as noted by the number of viable cells and their sub-types when only glycolytic components were varied. This, in part, explains why glycolytic bottlenecks in stem cell-derived -cells show non-robust glucose-stimulated-insulin-secretion behaviour [44, 45]. The calibrated PoMs generated in this study were able to capture electrical oscillations at low glucose, which is a unique feature of human -cells (compared to rodents). Thus, the methodology developed should be applicable for creating human-specifc 3D models of coupled -cell clusters with and without partner cells such as -, -, - and pancreatic polypeptide cells. When coupled with experimental studies, these 3D models could be invaluable to understand the efect of heterogeneity on intact islet function, and to develop sophisticated, robust tissue-engineered platforms for therapeutics.

Here we have demonstrated the use of experimental constraints to frst construct and then calibrate PoMs of three human cell types. The cardiomyocyte PoMs partially reproduced observable variation in the functional features for their respective cell types, albeit with clear systematic inaccuracies. The ability of the human -cell PoM to recapitulate the range of function observed in experimental -cell recordings is encouraging for applying this strategy to integrated function of -cells in human islets. Constructing new methods for constraining these PoM approaches will require continued advancement of single cell assays.

4 Functional Heterogeneity 61

### **References**

pluripotent stem cell derived cardiomyocytes: a dynamic clamp study with virtual IK1. *Frontiers in Physiology*, 6, 2015.

stem cell-derived β cells in vitro is inhibited by a bottleneck in glycolysis. *Cell Reports*, 31(6):107623, May 2020.

45. Ishan Goswami, Eleonora de Klerk, Phichitpol Carnese, Matthias Hebrok, and Kevin E Healy. Multiplexed microfuidic platform for stem-cell derived pancreatic islet cells. *Lab on a Chip*, 2022.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 5 Realizing Synaptic Signal Transmission During Astrocyte-Neuron Interactions within the EMI Framework**

Julia Gorman, Konstantin Holzhausen, Joyce Reimer, and Jørgen Riseth

**Abstract** The tripartite synapse or "neural threesome" refers to the interplay in the synapse between neighbouring neurons, the synaptic cleft, and the surrounding glial cells. Despite extensive research, the efects of glial cells, such as astrocytes, on signal transduction between neurons are not fully understood. The Kirchhof-Nernst-Planck (KNP) and Extracellular-Membrane-Intracellular (EMI) models constitute a promising framework for modeling these kinds of systems. However, they lack the neurotransmitter-related mechanisms that are necessary to bridge signal transduction across the synaptic cleft. Here, we propose an extension to the KNP-EMI model by a spatio-temporal difusion-based description of the most prominent neurotransmitter, glutamate, that allows for investigation of the contribution of astrocytes to the functionality of the synapse. We validate our model by showing that the presence of an astrocyte in the domain afects the glutamate fux across the postsynaptic terminal, as observed physiologically. The proposed extension ofers a sufciently simple way of integrating synaptic glutamate dynamics into the KNP-EMI framework. It introduces the relevant interactions between electrical activity and difusion processes at the tripartite synapse that are necessary to assess how astrocytes might contribute to the functionality of the synapse. This work has implications for future studies involving glial mechanisms and other charged species within the KNP-EMI framework.

Jørgen Riseth (corresponding author)

Julia Gorman Department of Neurosciences, University of California, San Diego, USA

Konstantin Holzhausen Department of Physics, University of Oslo, NO

Joyce Reimer Division of Biomedical Engineering, University of Saskatchewan, CA

Department of Numerical Analysis and Scientifc Computing, Simula Research Laboratory, NO e-mail<EMAIL_ADDRESS>
### **5.1 Introduction**

Circuits in the brain are composed of complex connections of both neurons and glia. Neurons communicate via patterns of electrical signals that are propagated by molecules and ions. These signals, called action potentials, communicate information that gives rise to cognitive function, behaviors, and movements. Action potentials travel from neuron to neuron across small gaps called synapses. The neuron sending the signal, the presynaptic neuron, releases specifc molecules called neurotransmitters to the receiving neuron, or the postsynaptic neuron. However, neurons are not the only cells that are found in the synapse. A type of glial cell called an astrocyte envelopes the synapse and alters information transmission by modulating the neurotransmitters as they difuse across the synapse [1]. Thus, together, the neurons, the synapse, and the astrocyte form a "neural threesome" [2], or tripartite synapse. The role of this astrocyte mediated modulation in electrical and chemical interplay is not well understood.

One of the ways astrocytes are thought to afect synaptic transmission is through reuptake of the neurotransmitter glutamate. Once an electrical signal reaches the end of the presynaptic neuron, prepackaged vesicles release glutamate into the synaptic cleft. After release, glutamate difuses across the synapse where it binds AMPA receptors on the postsynaptic neuron, resulting in a net inward sodium (Na+ ) current through the ligand-gated AMPA receptor ion channels. As a result, the intracellular potential starts to increase which triggers the opening of nearby voltage-gated Na+ channels. This channel opening allows for a larger Na+ infux, bringing the postsynaptic neuron closer to its threshold for depolarization. If the signal is strong enough, subsequent ion dynamics across the postsynaptic membranes result in the onset of a new action potential and signal propagation continues.

As glutamate difuses across the synapse, astrocytes take up a percentage of the glutamate ions through high-afnity glutamate transporters [3, 4]. This mechanism is suggested to be a way of reducing cellular toxicity, as too much glutamate in the synapse can have excitotoxic efects [5]. Recently, it has been found that astrocytes also release glutamate to neurons in response to extracellular glutamate cues [6]. Because of this dual mechanism of glutamate uptake and release, astrocytes act as a bufer to balance excitatory and inhibitory neuronal activity.

Here, we propose a simple mechanistic model of glutamate in the tripartite synapse that allows for the assessment of these questions within the Kirchhof-Nernst-Planck (KNP) - Extracellular-Membrane-Intracellular (EMI) modeling framework. The KNP-EMI model is a mechanistic mathematical model for excitable tissue in spatiotemporal resolution. It is stated in the form of a set of coupled partial diferential equations, and includes interactions between ionic transmembrane currents and electric potentials [7, 8, 9]. In [7], Ellingsrud *et al.* show that their KNP-EMI model accurately describes depolarization along excitable membranes. Thus, it constitutes a truly spatio-temporal model of signal transduction along a cell, for example, along an axon or a dendrite of a neuron. However, the model in its current form does not account for neurotransmitters in the extracellular space, nor for neurotransmitterrelated mechanisms and interactions. For this reason, the KNP-EMI model is not able to capture the regulation of signal transduction at the synapse.

We propose an extension by introducing glutamate to the model. In order to bridge the synaptic gap, we identify three relevant biophysical mechanisms: glutamate release into the cleft by exocytosis from the presynaptic neuron, difusion in the extracellular space, and binding to postsynaptic AMPA receptors. All three of these processes are considered in our implementation. In addition to this, we assume that glutamate uptake is the fundamental interaction of the astrocyte at the synapse whereby astrocytes regulate signal transmission. Therefore, our proposed glutamate difusion model describes the properties of signal transduction unique to the synapse; in particular, how astrocytes afect the glutamate gradient and how a changing glutamate gradient afects signal propagation. Our extension to the KNP-EMI model features a simplifed yet qualitatively complete description of signal transmission across neurons. The extended model will allow for a comprehensive investigation of the astrocyte's role in signal propagation across neurons in a way that takes spatial efects, such as explicit geometries and local concentrations, into account.

The EMI model constitutes a general framework for modeling excitable tissue ranging from cardiac cells to neural systems [10]. Besides its broad range of application, it allows for the incorporation of explicit morphologies within the system of interest. In combination with the KNP framework, the KNP-EMI model allows for the representation of a complicated geometry and detailed electrochemical profles [8]. In order to leverage its power, we design our glutamate model under consideration of its integrability into the KNP-EMI model. Although its formulation is independent of the respective morphology, our implementation assumes a simplistic geometry in order to show its functionality as a proof of concept.

Our paper is structured as follows. First, we familiarize the reader with the geometry of the problem conceptually and derive the glutamate model. Subsequently, we demonstrate a way to estimate the model parameters based on literature values, restricting our model to hippocampal synapses. After providing implementation details, we validate the glutamate model by showing that it behaves as intended. Finally, we discuss the implications of our model within the context of the KNP-EMI framework.

### **5.2 Methods**

The glutamate model consists of two parts. The frst part describes the transport of glutamate across the synaptic cleft by difusion, whereas the second part describes how AMPA receptors react in response to glutamate to trigger a Na+ current across the postsynaptic terminal. Here, we describe the computational domain and the governing equations for the glutamate model. We then describe the parameters that are used within the model, and fnally, the numerical methods used in our simulations.

### **5.2.1 Representation of the Computational Domain**

The tripartite synapse features a variety of membranes that behave diferently. The characteristic behavior of each membrane is determined by the density and type of functional components it exhibits, e.g. channels, pumps, and receptors. Figure 5.1 shows a schematic of our envisioned computational domain for the full KNP-EMI model. It resembles a simplifed and abstracted two-dimensional realization of a tripartite synapse setup. The domain is divided into two distinct subdomains: the

Fig. 5.1: Representation of the computational domain in our glutamate KNP-EMI model. Synthesis of computational domains, channel dynamics, and boundary conditions.

extracellular space, Ωe, and the intracellular space, Ωi . Ωi represents the interior of the postsynaptic neuron's dendrite. The interface boundary separating Ωe from Ωi consists of the postsynaptic terminal, Γpost, and extrasynaptic membrane of the dendrite, Γd. We assume that Γd is impermeable to glutamate. In contrast, Γpost is a region with an afnity for glutamate due to the presence of AMPA receptors, which are only located on the part of the postsynaptic membrane that faces the synapse [11]. The exterior boundary consists of three biophysically distinct types. Γpre denotes the membrane between the presynaptic neuron and the synapse. Γa represents the membrane on the astrocyte, which also has an afnity for glutamate due to the presence of glutamate transporters that clear it from the synapse [4]. The remaining part of the exterior boundary is open for glutamate to pass, since the neurotransmitter freely difuses through the extracellular space out of the synaptic region. We note that because this report only covers the glutamate model, our computational domain is for now restricted to the extracellular space, Ω.

### **5.2.2 Mathematical Modeling of Glutamate Dynamics**

The glutamate transport is modelled using a pure difusion equation: fnd the glutamate concentration = (,) such that

$$
\partial\_t \lg \mathbf{g} = D\_\mathbf{e} \Lambda \mathbf{g} \quad \text{on} \quad \mathfrak{Q}\_\mathbf{e} \times [0, T] \quad . \tag{5.1}
$$

Here, is the difusion coefcient for glutamate within the extracellular space. The glutamate release from the presynaptic terminal and absorption of glutamate on the various cell membranes are modelled using the following combination of von Neumann and Robin boundary conditions:

$$-D\_{\mathbf{e}} \nabla \mathbf{g} \cdot \mathbf{n}\_{\Gamma\_{\mathbf{d}}} = \qquad \qquad \mathbf{0} \qquad \qquad \text{on} \quad \Gamma\_{\mathbf{d}} \times [0, T] \tag{5.2}$$

$$-D\_{\mathbf{e}} \nabla \mathbf{g} \cdot \mathbf{n}\_{\Gamma\_{\mathbf{a}}} = \begin{array}{c} k\_{\text{update}} \ \mathbf{g} \vert\_{\Gamma\_{\mathbf{a}}} \qquad \text{on} \quad \Gamma\_{\mathbf{a}} \times \begin{bmatrix} 0,T \end{bmatrix} \tag{5.3}$$

$$-D\_{\rm e} \nabla \mathbf{g} \cdot \mathbf{n}\_{\Gamma\_{\rm post}} = \quad k\_{\rm bind} \ \mathbf{g} \vert\_{\Gamma\_{\rm post}} \qquad \text{on} \quad \Gamma\_{\rm post} \times \{0, T\} \tag{5.4}$$

$$-D\_{\rm e} \nabla g \cdot \mathbf{n}\_{\Gamma\_{\rm pre}} = -\sum\_{i} v\_{\nu}^{i} \,\delta\left(t - t\_{i}\right) \quad \text{on} \quad \Gamma\_{\rm pre} \times \{0, T\}. \tag{5.5}$$

Note that here, **n** denotes the outward pointing normal vectors of the respective boundary surfaces. As a consequence, the left hand sides denote the outward fuxes of glutamate through them. We model glutamate binding to the AMPA receptors as well as glutamate uptake at the astrocyte's membrane by frst order reaction kinetics. Thus, bind specifes the afnity of the AMPA receptor towards glutamate and uptake the absorption rate. The refecting boundary condition at Γd in Equation (5.2) realizes the impermeability of glutamate through the extrasynaptic membranes of the dendrite. Equation (5.5) describes the source of glutamate in our model, which is vesicular release due to exocytosis. In accordance with Clements *et al.* [12] and Jonas *et al.* [13], we model an instantaneous release of vesicles at release times . denotes the vesicular spatial profle defned in Equation (5.9). On the remaining open boundaries, we assume homogeneous Dirichlet boundary conditions

$$\|\mathbf{g}\|\_{\Gamma\_{\text{open}}} = \mathbf{0} \quad . \tag{5.6}$$

This choice is made for simplicity, and future work will beneft from more realistic modelling of free difusion across the open boundaries. Since glutamate is supposed to be cleared from the synaptic cleft after exocytosis, we assume an initial equilibrium of (, = 0) = 0.

### **5.2.3 Modeling AMPA Gating Dynamics**

We model the dynamics of the AMPA receptors on the postsynaptic terminal using a model proposed by Tewari and Majurmdar [14]. The ODE model describing the AMPA gating dynamics reads

$$\frac{d}{dt}\,m\_{\rm AMPA} = \alpha\_{\rm AMPA} \, g \Big|\_{\Gamma\_{\rm paa}} \left(1 - m\_{\rm AMPA}\right) - \beta\_{\rm AMPA} \, m\_{\rm AMPA} \quad . \tag{5.7}$$