myfnsymbols** ††

# Homogenisation for the monodomain model in the presence of microscopic
fibrotic structures

Brodie A. J. Lawson  Centre for Data Science, Queensland University of
Technology, Queensland, Australia  Rodrigo Weber dos Santos  Graduate Program
on Computational Modeling, Universidade de Juiz de Fora, Minas Gerais, Brazil
Ian W. Turner  ARC Centre of Excellence for Mathematical and Statistical
Frontiers, Queensland University of Technology, Queensland, Australia  Alfonso
Bueno-Orovio  Department of Computer Science, University of Oxford,
Oxfordshire, United Kingdom  Pamela Burrage  ARC Centre of Excellence for
Mathematical and Statistical Frontiers, Queensland University of Technology,
Queensland, Australia  Kevin Burrage  ARC Centre of Excellence for
Mathematical and Statistical Frontiers, Queensland University of Technology,
Queensland, Australia

###### Abstract

Computational models in cardiac electrophysiology are notorious for long
runtimes, restricting the numbers of nodes and mesh elements in the numerical
discretisations used for their solution. This makes it particularly
challenging to incorporate structural heterogeneities on small spatial scales,
preventing a full understanding of the critical arrhythmogenic effects of
conditions such as cardiac fibrosis. In this work, we explore the technique of
homogenisation by volume averaging for the inclusion of non-conductive micro-
structures into larger-scale cardiac meshes with minor computational overhead.
Importantly, our approach is not restricted to periodic patterns, enabling
homogenised models to represent, for example, the intricate patterns of
collagen deposition present in different types of fibrosis. We first highlight
the importance of appropriate boundary condition choice for the closure
problems that define the parameters of homogenised models. Then, we
demonstrate the technique’s ability to correctly upscale the effects of
fibrotic patterns with a spatial resolution of
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ into much larger numerical mesh
sizes of 100-$250\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The homogenised
models using these coarser meshes correctly predict critical pro-arrhythmic
effects of fibrosis, including slowed conduction, source/sink mismatch, and
stabilisation of re-entrant activation patterns. As such, this approach to
homogenisation represents a significant step towards whole organ simulations
that unravel the effects of microscopic cardiac tissue heterogeneities.

##  1 Introduction

Computational simulation plays a critical role in our understanding of the
functioning of the heart, in particular the complex manifestations of its
excitable media dynamics into dangerous arrhythmias [1]. An important
contributor to many types of arrhythmia is cardiac fibrosis, the pathological
formation of scar tissue in the heart [2], its arrhythmogenic impacts
depending on its spatial organisation on both microscopic [3] and macroscopic
[4] scales. However, owing both to limitations of computational feasibility
and the resolution of clinical imaging approaches, anatomically-accurate
meshes used for the simulation of electrical signalling in the heart typically
have spacings of minimum 100 micrometres [5]. This is at least an order of
magnitude too large to resolve the complex and varied microscopic structures
of fibroblast-deposited collagen that interfere with wave propagation [6]. It
is therefore vital that such “sub-mesh scale” effects of fibrotic obstacles be
incorporated into simulations without altering the mesh spacing. Even when not
working with an anatomical mesh, this type of upscaling represents a
significant computational time saving that may be used alongside other
acceleration techniques such as improved numerical techniques (e.g. [7])
and/or hardware architectures [8].

There has been some progress in incorporating small-scale fibrotic structures
into larger-scale cardiac electrophysiology simulations. Through a clever node
re-labelling, Costa et al. were able to incorporate disconnections between
neighbouring elements due to strands of collagen [9], although such an
approach does not necessarily account for the effects of obstacles within mesh
elements. An alternative approach is based on the mathematical technique of
homogenisation, which explicitly seeks to represent micro-scale effects as
modifications to macro-scale problems [10]. This has some history in cardiac
electrophysiology in the derivation of the well-known bidomain model [11, 12,
13, 14, 15], or its modification in the case of less ordered arrangements of
cells [16] or to represent the non-ohmic nature of tissue conduction [17].
However, research in using homogenisation to incorporate the effects of
fibrosis have been largely limited to spatially periodic structures [18, 19].
Austin et al. [20] used homogenisation based on multigrid techniques to
incorporate arbitrarily arranged obstacles into larger-scale simulations, but
the analysis of spatiotemporal dynamics was not given. Most importantly, none
of these approaches have considered whether or not homogenisation is able to
capture the mechanisms through which microscopic obstacles to conduction act
as arrhythmia precursors. This is arguably the primary aspect of interest in
cardiac electrophysiology simulations.

In this work, we use a volume averaging approach for the incorporation of
arbitrary structures of microscopic obstacles into a larger-scale problem. We
explore several different choices of boundary conditions for homogenisation
sub-problems, in order to determine which is most appropriate for the
challenging case of sharp-fronted travelling wave dynamics in the presence of
completely non-conductive obstacles of arbitrary shapes. We demonstrate the
successful capture of several important pro-arrhythmic effects of cardiac
fibrosis by block homogenised models, with one to three orders of magnitude
fewer nodes than the corresponding fine-scale models. Indeed, some type of
homogenisation is likely necessary for the inclusion of small-scale fibrotic
structures into a typical three-dimensional mesh of even a single heart
chamber, owing to the number of nodes/elements that would compose the
corresponding fine-scale discretisation.

## Materials and Methods

### Simulation of cardiac excitation in obstructed tissue

The dynamics of cardiac excitation are here governed by the monodomain model
[21], a simplification of the bidomain model that offers similar quality of
predictions in many contexts [22, 23, 24]. The monodomain model is a parabolic
partial differential equation coupled to a set of ordinary differential
equations. In the presence of non-conductive obstacles, the monodomain model
may be expressed in the form

$\displaystyle\frac{\partial v}{\partial t}$
$\displaystyle=\nabla\cdot\Bigl{(}\mathbf{D}\nabla
v\Bigr{)}-\frac{1}{C_{m}}(I_{\mbox{\scriptsize
ion}}(v,\mathbf{s})+I_{\mbox{\scriptsize stim}})$ within conducting tissue (1)
$\displaystyle\frac{d\mathbf{s}}{dt}$ $\displaystyle=\mathbf{f}(v,\mathbf{s})$
within conducting tissue $\displaystyle 0$ $\displaystyle=(\mathbf{D}\nabla
v)\cdot\hat{\mathbf{n}}$ $\displaystyle\mbox{on boundaries (incl.
obstructions)}.$

Here $v$ is the membrane potential (in lower case to denote a micro-scale
variable), $C_{m}$ is the membrane capacitance and $\mathbf{D}$ is the
conductivity tensor. $I_{\mbox{\scriptsize stim}}$ refers to externally
supplied stimulus current, and $I_{\mbox{\scriptsize ion}}$ specifies the flow
of ions in/out of cardiac cells, which depends on both the membrane potential
and a set of state variables $\mathbf{s}$. We choose the reduced version of
the ten Tusscher et al. ionic model [25] to define $I_{\mbox{\scriptsize
ion}}$ and $f$. This model represents action potentials in human ventricular
epicardium using formulations for all of the major Na+, K+, Ca2+ currents
involved, but with a set of simplifying assumptions made to greatly reduce
computational cost.

Fibrotic obstructions are defined on a fine-scale grid of spacing $\Delta
x=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, a similar order to the pixels
in histological images indicating the spatial arrangement of collagenous
obstacles in cardiac fibrosis [6]. The homogenised models we construct seek to
represent the effects of these obstacles on a regular grid of a larger scale,
as visualised in Fig 1. We consider homogenised models with $\Delta
x=50,100,250,$500\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, corresponding to
grid spacings common in computational cardiac electrophysiology. Although the
current setup uses regular grids with the edges of large scale mesh elements
aligned with finescale grid elements, for simplicity, the homogenisation
theory presented subsequently does not depend on these choices.

Figure 1: The mesh and sub-mesh involved with homogenisation. An example
pattern of collagenous obstructions due to fibrosis (dark red) defined on a
small-scale mesh, and the larger mesh on which effective conductivity tensors
are defined. A single averaging volume ($\Omega$) is marked in black, along
with its associated conductive region, $\Omega_{m}$, in orange. The black
dotted line shows the region on which the closure problem associated with
$\Omega$ is solved when a layer of skin is included around the averaging
volume (see Methods and Materials). An example control volume used by the
finite volume approach used for numerical discretisation is also pictured, in
blue.

Eq (1), the homogenised equivalent and the closure problems described
subsequently, are all solved using a vertex-centred finite volume method.
Integrals over control volumes are approximated by exactly integrating linear
interpolants, constructed using the nodal values at element vertices. This
choice results in a numerical scheme that is similar to that derived using a
finite element method without mass lumping, an approach that has been shown to
reduce the sensitivity of the monodomain model to the grid spacing used [5].
Timestepping is performed using the second-order generalisation of the
Rush–Larsen method presented by Perego and Veneziani [26], with a timestep of
$\Delta t=$0.05\text{\,}\mathrm{ms}$$.

Minimisation of sensitivity of the monodomain model to the spatial
discretisation used is paramount for both the evaluation and utilisation of
homogenisation, as converting to a homogenised large-scale problem of course
incurs a significant change in node spacing. As such, we further correct for
the effect of the grid spacing by multiplying all conductivity tensors in our
homogenised problems by a constant, such that the conduction velocity in a
one-dimensional (non-fibrotic) fibre is consistent with that predicted using
the finescale grid spacing ($10\text{\,}\mathrm{\SIUnitSymbolMicro m}$). We
note that this correction factor is selected before, and entirely independent
from, the homogenisation process and thus does not act to inflate the
perceived performance of the homogenisation itself.

### Homogenised monodomain models

Homogenisation is approached via the method of volume averaging, detailed in
[27]. Volume averaging, in the case where non-conductive material is present,
makes use of a pair of averaging operators that average a quantity over the
conductive portions, $\Omega_{m}$, of an averaging volume, $\Omega$. Here, as
is typical for this kind of homogenisation, the averaging volume is taken as a
single element of the large-scale grid (see Fig 1). The two averaging
operators are

$\displaystyle\langle\cdot\rangle=\frac{1}{|\Omega_{m}|}\int\limits_{\Omega_{m}}\cdot\,d\Omega_{m},$
intrinsic average $\displaystyle\langle\cdot\rangle_{\mbox{\tiny
sup}}=\frac{1}{|\Omega|}\int\limits_{\Omega_{m}}\cdot\,d\Omega_{m},$
$\displaystyle\mbox{superficial average},$

and are thus linked by the volume fraction of conductive material,
$\phi=|\Omega_{m}|/|\Omega|$, as

$\langle\cdot\rangle_{\mbox{\tiny sup}}=\phi\langle\cdot\rangle.$

We present a brief, simplified derivation of how these operators are used to
formulate homogenised versions of Eq (1).

Consider, as an example, an averaging volume occupied by only a small
proportion of cardiac tissue with the remainder occupied by fibrotic
obstruction. Even if all of the conductive material on the finescale has a
membrane potential $v$ that is above the depolarisation threshold and is about
to excite, $\langle v\rangle_{\mbox{\tiny sup}}$ could very well fall below
this threshold and hence fails to describe even the “average” dynamics of the
tissue in the averaging volume. As such, the intrinsic average $\langle
v\rangle$ is the quantity we must use in our homogenised model. However, the
superficial averaging operator allows us to use the spatial averaging theorem
[27, Ch. 1],

$\biggl{\langle}\nabla\cdot\mathbf{f}\biggr{\rangle}_{\mbox{\tiny
sup}}=\nabla\cdot\biggl{\langle}\mathbf{f}\biggr{\rangle}_{\mbox{\tiny
sup}}+\frac{1}{|\Omega|}\int_{\sigma_{b}}\mathbf{f}\cdot\hat{\mathbf{n}}\,d\sigma_{b}$
(2)

where $\sigma_{b}$ denotes the surface of the interface between the conductive
and non-conducting regions.

Applying the superficial averaging operator to both sides of the spatially-
varying equation in Eq (1), and defining $J=(I_{\mbox{\scriptsize
ion}}(v,\mathbf{s})+I_{\mbox{\scriptsize stim}})/C_{m}$ to simplify notation,
then

$\left\langle\frac{\partial v}{\partial t}\right\rangle_{\mbox{\tiny
sup}}=\biggl{\langle}\nabla\cdot\Bigl{(}\mathbf{D}\nabla
v\Bigr{)}\biggr{\rangle}_{\mbox{\tiny
sup}}-\Bigl{\langle}J(v,\mathbf{s})\Bigr{\rangle}_{\mbox{\tiny sup}}.$ (3)

Using Eq (2) gives

$\left\langle\frac{\partial v}{\partial t}\right\rangle_{\mbox{\tiny
sup}}=\nabla\cdot\Bigl{\langle}\mathbf{D}\nabla
v\Bigr{\rangle}_{\mbox{\scriptsize
sup}}+\frac{1}{|\Omega|}\int\limits_{\sigma_{b}}(\mathbf{D}\nabla
v)\cdot\hat{\mathbf{n}}\,d\sigma_{b}+\Bigl{\langle}J(v,\mathbf{s})\Bigr{\rangle}_{\mbox{\tiny
sup}}$

however the integral term is seen to be zero due to the boundary condition for
the interface between conductive tissue and obstacles. Additionally, in the
absence of considering electromechanical coupling we have that obstacles are
fixed and so the order of averaging and time differentiation on the left hand
side may be reversed. Together, these observations give

$\frac{\partial\langle v\rangle_{\mbox{\tiny sup}}}{\partial
t}=\nabla\cdot\biggl{(}\Bigl{\langle}\mathbf{D}\nabla
v\Bigr{\rangle}_{\mbox{\tiny
sup}}\biggr{)}+\Bigl{\langle}J(v,\mathbf{s})\Bigr{\rangle}_{\mbox{\tiny sup}}$

or in terms of the intrinsic average,

$\phi\frac{\partial\langle v\rangle}{\partial
t}=\nabla\cdot\biggl{(}\phi\Bigl{\langle}\mathbf{D}\nabla
v\Bigr{\rangle}\biggr{)}+\phi\Bigl{\langle}J(v,\mathbf{s})\Bigr{\rangle}.$ (4)

Finally, we wish to express Eq (4) solely in terms of a macroscopic variable,
$V=\langle v\rangle$. In order to achieve this, we define the effective
conductivity tensor such that

$\mathbf{D}_{\mbox{\tiny eff}}\nabla V=\Bigl{\langle}\mathbf{D}\nabla
v\Bigr{\rangle}$ (5)

and make the simplifying approximation

$\Bigl{\langle}J(v,\mathbf{s})\Bigr{\rangle}\approx J(V,\mathbf{s}).$ (6)

This results in the homogenised monodomain model describing the large-scale
behaviour of the system,

$\displaystyle\phi\frac{\partial V}{\partial t}$
$\displaystyle=\nabla\cdot\biggl{(}\phi\mathbf{D}_{\mbox{\tiny eff}}\nabla
V\biggr{)}+\phi J(V,\mathbf{s})$ within tissue (7)
$\displaystyle\frac{d\mathbf{s}}{dt}$ $\displaystyle=\mathbf{f}(V,\mathbf{s})$
within tissue $\displaystyle 0$ $\displaystyle=(\mathbf{D}_{\mbox{\tiny
eff}}\nabla V)\cdot\hat{\mathbf{n}}$ $\displaystyle\mbox{on problem
boundaries}.$

Note that obstacles no longer act through the boundary conditions, but instead
through their effect on $\mathbf{D}_{\mbox{\tiny eff}}$ and $\phi$.

We briefly note that a separate technique, the smoothed boundary method, also
shifts such boundary conditions into the governing equation to arrive at the
formulation Eq (7) [28, 29]. Volume averaging theory thus serves as a robust
means of deriving the smoothed boundary approach. The key difference between
the two approaches is their use cases. Smoothed boundary methods use
especially fine grids at boundaries in order to accurately resolve their
effects, whereas homogenisation by volume averaging instead seeks to represent
these effects only on a larger scale, typically to greatly reduce
computational demand.

### Determination of effective conductivities

Eq (5) relates the macroscopic gradient of potential $\nabla V$ to the fine-
scale gradient, $\nabla v$. These gradients will change through the course of
a simulation of Eq (1) or Eq (7), but fixed conduction tensors can be set by
instead solving separate subproblems in which a macroscopic gradient is
artificially applied [30]. Imposing a macroscopic unit gradient in the $i$-th
direction ($\nabla V=\mathbf{e}_{i}$) over the averaging volume, the $i$-th
column of the effective conductivity tensor is then given by

$\mathbf{D}_{\mbox{\tiny eff}}\mathbf{e}_{i}=\Bigl{\langle}\mathbf{D}\nabla
v_{i}\Bigr{\rangle},\qquad\qquad i=1,...,d,$ (8)

with $\mathbf{e}_{i}$ the standard basis vectors in the $d$-dimensional space.
Eq (8) applies regardless of the shape of the averaging volume used, or
phrased differently, any set of applied gradients may be used to calculate the
elements of $\mathbf{D}_{\mbox{\tiny eff}}$ so long as they are in linearly
independent directions [30].

Each $v_{i}$ is the solution of a closure problem, a micro-scale subproblem
determining how the imposed macroscopic gradient translates into flow through
the material being homogenised (according to the diffusive portion of Eq (1)).
These problems (and their boundary conditions discussed subsequently) are more
naturally expressed in terms of corresponding closure variables,

$w_{i}=v_{i}-x_{i},$ (9)

the substitution of which into Eq (1) results in microscale transport being
defined by

$\displaystyle 0$ $\displaystyle=\nabla\cdot\Bigl{(}\mathbf{D}(\nabla
w_{i}+\mathbf{e}_{i})\Bigr{)}$ within conductive tissue (10) $\displaystyle 0$
$\displaystyle=\Bigl{(}\mathbf{D}(\nabla
w_{i}+\mathbf{e}_{i})\Bigr{)}\cdot\hat{\mathbf{n}}$ $\displaystyle\mbox{on
boundaries with obstructions}.$

Solutions of Eq (10) for each choice of $\mathbf{e}_{i}$ define the different
$w_{i}$’s, and hence $v_{i}$’s, with which each column of a large scale mesh
element’s effective conductivity tensor is then calculated by Eq (8). Each
element in the large scale problem is assigned its own conductivity tensor via
the solution of its own set of independent closure problems, a method referred
to as block homogenisation [31].

### Closure problem boundary conditions

Arbitrary structures of obstructions in cardiac tissue, such as collagen in
cardiac fibrosis, vary significantly in terms of both patterning [6] and
density (for example, within or outside of an infarct region). As such, the
assumption of periodicity used in many homogenisation approaches may not be
inappropriate. In particular, this assumption will not function when
conduction takes place through thin channels stretching across different
averaging volumes, as unless those channels happen to align at the opposite
ends of an individual volume, the periodic extension implies a non-conducting
structure [10]. This scenario arises in the cardiac electrophysiology context
through conductive isthmuses (channels) running through non-conductive scar
regions, which are of particular interest as potential substrates for
arrhythmia [32].

To address this, we consider alternative boundary conditions for closure
problems, Eq (10), that relax the assumption of periodicity and are commonly
used in Laplacian homogenisation [30, 33]. For a rectangular averaging volume
with side lengths $L_{i}$, the set of boundary conditions considered are
expressed

$\displaystyle Periodic:$ $\displaystyle
w_{i}(\mathbf{x}+L_{j}\mathbf{e}_{j})$ $\displaystyle=w_{i}(\mathbf{x})$
$\displaystyle\forall j\in{1,\ldots,d}$ (11) $\displaystyle Linear:$
$\displaystyle w_{i}$ $\displaystyle=0$
$\displaystyle\mathbf{x}\in\partial\Omega$ $\displaystyle Confined:$
$\displaystyle w_{i}$ $\displaystyle=0$
$\displaystyle\mathbf{x}\in\partial\Omega,\>\mathbf{e}_{i}\cdot\hat{\mathbf{n}}\neq
0$ $\displaystyle(\nabla w_{i}+\mathbf{e}_{i})\cdot\hat{\mathbf{n}}$
$\displaystyle=0$
$\displaystyle\mathbf{x}\in\partial\Omega,\>\mathbf{e}_{i}\cdot\hat{\mathbf{n}}=0.$

Recalling that the closure problems operate by imposing a macroscopic gradient
and considering the resulting flow through the averaging volume, linear
boundary conditions may be interpreted as holding all boundaries fixed
according to the imposed gradient. Confined conditions, where two opposing
boundaries are fixed to maintain the imposed gradient while the remaining
boundaries are given no-flux conditions, make Eq (8-10) a numerical recreation
of Darcy’s experiments that first derived hydraulic conductivity [34].

To further reduce the effects of the assumptions implied by different choices
of boundary conditions, we also consider the effects of including a layer of
“skin” around the averaging volume. This approach extends the domain on which
Eq (10) is solved for each averaging volume, beyond the boundaries of that
averaging volume [35, 36] (see also Fig 1). As the calculation of effective
conductivity through Eq (8) still only considers the averaging volume, this
has the effect of moving the boundaries away from the region used for
calculation and thus hopefully minimising boundary condition effects. This is
strongly related to the concept of “oversampling” in other homogenisation
contexts [37].

The downside of including skin around averaging volumes is the loss of the
guarantee that effective conductivity tensors are symmetric for linear or
periodic boundary conditions [30]. To account for this (as well as asymmetric
tensors that may be produced when using confined boundary conditions), we use
the algorithm of Higham [38] to find the symmetric semipositive definite
tensor that is closest (in terms of Frobenius norm) to the calculated
effective conductivity, $\mathbf{D}_{\mbox{\tiny eff}}$. Briefly, this
approach works by calculating the eigendecomposition of the symmetric portion
of the initial tensor, $(\mathbf{D}_{\mbox{\tiny eff}}+\mathbf{D}_{\mbox{\tiny
eff}}^{T})/2$, zeroing out its negative eigenvalues, and then rebuilding it
using this modified eigendecomposition.

## Results

### Non-conductive obstacles complicate boundary condition selection

The performance of different boundary conditions in homogenisation is often
discussed in reviews of the method [34, 33]. However, the case where wholly
non-conductive material is present and extends to a significant portion of the
boundaries appears to be far less well-considered. We address this by
exploring how the three different choices of boundary conditions perform on
two particularly illustrative examples, highlighting just how critical the
choice of boundary conditions for homogenisation subproblems can be. This then
also informs the choice of boundary conditions for formulating homogenised
monodomain models in different contexts, as well as the interpretation of the
results we observe for the homogenised monodomain models we consider here.

#### Thin Barriers

The first scenario we consider is a single large-scale element composed of
isotropically conductive material of constant property, but with a thin strip
of obstructive material running vertically along its whole length. This
structure is pictured in Fig 2, and it can be intuited that such a structure
permits no macroscopic flow in the $x$-direction. In the $y$-direction,
although the amount of macroscopic transport will be lowered slightly by the
non-conductive portion, the structure poses no obstacle to vertical
conduction. As such, our “intrinsic” formulation should define the effective
conductivity for this element to be

$\mathbf{D}_{\mbox{\tiny eff (true)}}=\begin{pmatrix}0&0\\\ 0&D\end{pmatrix},$

where $D$ is the (scalar) conductivity of the conductive medium. However, not
all choices of boundary conditions are able to obtain this simple result.

True | Linear BCs | Confined BCs | Periodic Conditions
---|---|---|---
$\mathbf{D}_{\mbox{\tiny eff}}=\begin{pmatrix}0&0\\\ 0&1\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny eff}}=\begin{pmatrix}0.486&0.000\\\ 0.000&1.000\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny eff}}=\begin{pmatrix}0.000&0.000\\\ 0.000&1.000\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny eff}}=\begin{pmatrix}0.000&0.000\\\ 0.000&1.000\end{pmatrix}$
| $\mathbf{D}_{\mbox{\tiny skin}}=\begin{pmatrix}0.253&0.000\\\ 0.000&1.000\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny skin}}=\begin{pmatrix}0.000&0.000\\\ 0.000&1.000\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny skin}}=\begin{pmatrix}0.000&0.000\\\ 0.000&1.000\end{pmatrix}$

Figure 2: Performance of different homogenisation boundary conditions (BCs)
for a pernicious test case. Conductive tissue is shown in yellow, and fibrotic
obstruction in dark red (dark grey in visualisations of the closure problem
solution). The dashed blue rectangle indicates the region averaged over when
considering the effects of skin. Pictured are the solutions to the closure
subproblem Eq (10) for the variable $w_{1}$, which when averaged together with
the $w_{2}$ solutions (not pictured) result in the conductivity tensors given
in the table. Confined and periodic conditions result in a constant gradient
solution, such that $\bigl{\langle}\mathbf{D}(\mathbf{e}_{1}+\nabla
w_{1})\bigr{\rangle}=\mathbf{0}$ and zero horizontal flow is correctly
predicted. Linear boundary conditions result in a degradation of that solution
near the boundaries, and this results in an effective tensor that still
permits some horizontal flow. Use of skin to shift averaging away from the
boundaries helps reduce this effect but falls far short of eliminating it.

In Fig 2 we present $\mathbf{D}_{\mbox{\tiny eff}}$ as predicted by the
different choices of boundary conditions, with and without including a layer
of skin. Also presented are the corresponding solutions of Eq (10) when the
macroscopic gradient is imposed in the $x$ direction ($w_{1}$). Confined and
periodic conditions produce $w_{1}$ solutions with a constant gradient from
left to right, and this gradient has the appropriate magnitude such that the
quantity $\Bigl{\langle}\mathbf{D}\nabla
v_{1}\Bigr{\rangle}=\Bigl{\langle}\mathbf{D}(\nabla
w_{1}+\mathbf{e}_{1})\Bigr{\rangle}=\mathbf{0}$. Hence, Eq (8) results in the
correct conductivity tensor, with zero macroscopic flow in the horizontal
direction.

On the other hand, when linear boundary conditions fix $w_{1}=0$ along all
boundaries, the constant gradient solution is lost (Fig 2). This results in an
incorrect effective conductivity tensor that permits considerable flow in the
horizontal direction. By avoiding averaging over the boundaries where the
solution is most degraded, including skin has a significant positive effect on
the conductivity tensor calculated using linear boundary conditions. Even
still, a considerable amount of horizontal flow is permitted, and the effect
of the vertical barrier on conduction is essentially lost.

#### Diagonally-oriented Channels

The second example we consider is a structure for which the scenario is
essentially completely reversed. This structure is a pair of diagonally
conducting channels through an otherwise non-conductive medium (Fig 3).
Diagonal structures resulting in anisotropy are a popular choice for
indicating the effects of the choice of homogenisation boundary conditions on
the effective tensors calculated (for example [39]), but scenarios considered
typically treat both materials as conductive. For the situation of conducting
channels through a non-conductive medium that we consider here, the effects of
boundary condition choice are even more pronounced.

We again treat the conductive material as isotropic and with constant
diffusivity $D$. Flow is unimpeded in the direction of the channel, and zero
in the direction perpendicular to the channel, and thus by rotational
arguments,

$\mathbf{D}_{\mbox{\tiny eff (true)}}=\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}D&0\\\
0&0\end{pmatrix}\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}=D\begin{pmatrix}\cos^{2}\theta&\cos\theta\sin\theta\\\
\cos\theta\sin\theta&\sin^{2}\theta\end{pmatrix}.$

Again, as we take an intrinsic formulation, the width of the channel does not
appear in this result. The macroscopic amount of transport is controlled by
the volume fraction $\phi$, while $\mathbf{D}_{\mbox{\tiny eff}}$ as we define
it describes the character of this transport.

True | Linear BCs | Confined BCs | Periodic Conditions
---|---|---|---
$\mathbf{D}_{\mbox{\tiny eff}}=\begin{pmatrix}0.5&0.5\\\ 0.5&0.5\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny eff}}=\begin{pmatrix}0.518&0.481\\\ 0.481&0.518\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny eff}}=\begin{pmatrix}0.000&0.000\\\ 0.000&0.000\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny eff}}=\begin{pmatrix}0.000&0.000\\\ 0.000&0.000\end{pmatrix}$
| $\mathbf{D}_{\mbox{\tiny skin}}=\begin{pmatrix}0.500&0.500\\\ 0.500&0.500\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny skin}}=\begin{pmatrix}0.000&0.000\\\ 0.000&0.000\end{pmatrix}$ | $\mathbf{D}_{\mbox{\tiny skin}}=\begin{pmatrix}0.000&0.000\\\ 0.000&0.000\end{pmatrix}$

Figure 3: Performance of different homogenisation boundary conditions (BCs)
for a second pernicious test case. Conductive tissue is shown in yellow, and
fibrotic obstruction in dark red (dark grey in the visualisation of the
closure problem solution). The dashed blue rectangle indicates the region
averaged over when considering the effects of skin. Pictured is the solution
to the closure problem Eq (10) for the variable $w_{1}$, and schematic
diagrams that indicate the failure of the other types of boundary conditions
in this scenario. Again, the evident patterning in the closure problem
solution with linear boundary conditions is disrupted at the boundaries.
However, disruption is minimal (particularly when skin is included) and the
calculated tensor approximates the true tensor. With confined boundary
conditions, no flow can pass from left to right when the top and bottom
boundaries are blocked (and analogously for when a vertical gradient is
imposed), and the calculated tensor is the zero tensor. For periodic
conditions, a zero tensor is obtained because the periodic extension of the
pattern is seen to be a non-conducting structure.

Fig 3 shows how the different choices of boundary conditions perform in this
scenario (with $\theta=\pi/4$ and $D=1$). The closure problem solution $w_{1}$
again shows disruption at the boundaries when using linear boundary
conditions, but the effect is minor in this case and the effective tensor
successfully approximates the true effective conductivity. Using skin further
reduces the disruptive effect of the boundary condition, improving the
accuracy of the calculated tensor from one decimal place to three.

In contrast, periodic and confined boundary conditions now incorrectly predict
zero transport through the averaging volume, with the solution to Eq (10)
giving $\Bigl{\langle}\mathbf{D}(\nabla
w_{i}+\mathbf{e}_{i})\Bigr{\rangle}=\mathbf{0}$. This result is perhaps better
understood by considering the physical interpretation of this closure
equation, as presented in Fig 3. Confined boundary conditions imply a Darcy
experiment, and thus zero conductivity is predicted as no path leads through
the element from left to right or bottom to top. Periodic conditions result in
the channel being blocked at both ends, as demonstrated by the periodic
extension pictured in the figure. In the case where the diagonal channels are
positioned so that their beginnings and ends align in the periodic extension,
periodic conditions can then predict the correct tensor (results not
pictured).

Together with the results from the previous section regarding thin barriers,
we see that when the arrangement of non-conducting obstacles is arbitrary,
there exist scenarios for which each type of boundary conditions we have
considered produces a poor estimate of macroscopic conductivity. As such, we
cannot select a consistently superior choice and instead now consider their
performance for practical use of homogenisation in the context of cardiac
electrophysiology. Specifically, we explore the potential of homogenised
models to capture several key pro-arrhythmic effects of cardiac fibrosis, as
represented by the presence of non-conductive fibrotic deposits.

### Homogenised monodomain models capture macroscopic excitation propagation
in obstructed tissue

Fibrotic obstructions slow the propagation of cardiac excitation through
afflicted tissue, a key component of fibrosis’ pro-arrhythmic effect [40] as
it decreases the “wavelength” that governs the survival of dangerous re-
entries [41]. We therefore use the wavespeed through obstructed tissue as the
first test of our homogenisation approach, specifically two-dimensional slices
of cardiac tissue measuring $5\text{\,}\mathrm{cm}$ $\times$
$0.5\text{\,}\mathrm{cm}$ with an anisotropic conductivity tensor with faster
conduction (3:1 ratio for conductivity between the $x$ and $y$ directions) in
the direction of propagation to match the faster conduction along cardiac
fibres in the heart. In these slices of tissue, we place non-conductive
obstructions at random, either $10\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$\,\\!\times\,\\!$$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ or
$90\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$\,\\!\times\,\\!$$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in size, with
the latter oriented both parallel and perpendicular to the direction of
propagation (depicted in Fig 4). The effect of these three types of fibrosis
on conduction has also been recently considered, separate from the context of
homogenisation [42].

Figure 4:  Prediction by homogenised models of conduction speeds in different
types of obstructed tissue. Top Row: Representative sections of the
$5\text{\,}\mathrm{cm}$ $\times$ $0.5\text{\,}\mathrm{cm}$ fibres showing the
patterns of obstruction considered (obstacles in red). Pictured examples are
the case of 25% obstacles. Waves of excitation move from left to right. Bottom
Rows: Performance of the different types of boundary conditions on the three
types of obstruction, for different choices of averaging volume size and
boundary conditions for closure subproblems Eq (10). Where a dot does not
appear for a given level of obstruction, this corresponds to a failure to
propagate the length of the fibre. Both boundary condition selection and
averaging volume size have an important effect on homogenisation performance,
with periodic conditions most accurate for smaller averaging volumes $(\Delta
x\leq$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$)$ and linear conditions most
accurate for the larger ($\Delta x=$250\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$). Homogenised models with $\Delta
x=$500\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ are universally poor for the
more challenging, highly obstructed problems. Best overall performance is
obtained by using linear boundary conditions and a 25$\\!\times\\!$25
averaging volume ($10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ up to
$250\text{\,}\mathrm{\SIUnitSymbolMicro m}$). Particularly notable is the case
of “perpendicular” obstructions, where homogenised models using boundary
conditions other than linear are prone to over-predicting conduction block.

Wavespeeds predicted by the homogenised models match well with fine-scale
wavespeeds, but begin to deviate as the amount of fibrotic obstruction
increases and paths of conduction become more torturous (Fig 4). The
performance of different choices of boundary conditions for closure
subproblems is comparable, with the superior choice also depending on the size
of the averaging volume used. This is a result of the interaction between the
error due to the homogenisation, and the numerical consequences of changing
the grid spacing (which is exacerbated by the smaller conductivity tensors in
highly-obstructed tissue). Error is consistently worst for the largest
averaging volumes ($\Delta x=$500\text{\,}\mathrm{\SIUnitSymbolMicro m}$$),
the case in which homogenisation error is expected to be lowest as the ratio
between characteristic length scales grows smaller [27]. As such, this implies
that the effect of changing the grid spacing is the predominant source of
error.

Perpendicularly aligned fibrosis presents a particularly interesting scenario,
as it results in significant reduction in velocity even for small amounts of
obstruction, eventually culminating in complete block of conduction when the
proportion of obstructive material reaches 50%. Homogenised models will only
be able to predict this block if they feature conductivities small enough to
halt conduction, as all path information is lost. The results for this
perpendicular fibrosis case mirror those seen in the diagonal channel example
we consider in detail (Fig 3). Specifically, homogenised models attained using
linear boundary conditions can potentially over-predict successful conduction,
as a connection between any averaging volume boundaries will result in a
weakly conductive element (Fig 2) even when the fine-scale structure in fact
creates a dead end. On the other hand, confined and periodic boundary
conditions result in homogenised models that significantly over-predict
conduction block as the homogenisation process can only “see” conductive paths
that fit with the assumptions that underlie them. Overall, best performance is
seen using linear boundary conditions and $\Delta
x=$250\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, a choice that performs very
well across all of the different patternings of fibrosis tested.

In order to further explore the interaction between homogenisation and
numerical error, we also separately consider the results of fine-scale models
that have their conductivity fields $\mathbf{D}(\mathbf{x})$ replaced with the
effective fields $\mathbf{D}_{\mbox{\tiny eff}}(\mathbf{x})$ and volume
fractions $\phi(\mathbf{x})$ obtained through the different block
homogenisation approaches. This removes the numerical effects of changing the
grid spacing and allows homogenisation error to be more directly examined, but
of course does not represent a practical use of homogenisation as there is no
computational saving. The results of these tests (Fig 5) present two important
conclusions. Firstly, homogenised models perform very well overall, confirming
that it is their changed gridsize that produces most of their discrepancy from
the equivalent finescale model. Secondly, the overall best-performing boundary
conditions now switch to periodic boundary conditions. This is not so
surprising, as the tissue fibres here have a consistent patterning of
obstacles throughout, and it is therefore reasonable to take a representative
volume element and treat the medium as periodic [27]. This shows that in the
context of the monodomain model (or other models sensitive to the spatial
discretisation used to simulate them), the potentially compensatory balance
between homogenisation and grid error must be considered in selecting the
parameters of the homogenised model.

Figure 5:  Performance of homogenisation for propagation through obstructed
tissue when separated from numerical effects. Displayed conduction velocities
are those predicted by homogenised models Eq (7) with different choices of
boundary conditions (BCs) and averaging volume size, but solved on the
original finescale grid to remove the effects of changing the numerical
discretisation. All homogenised models perform well, with best performance
obtained using the largest averaging volumes to maximise the difference in
length scales between finescale features and the homogenised model. Periodic
BCs perform best overall, as the medium is of consistent property and hence
admits a representative volume element. When obstacles are aligned
perpendicular to the direction of propagation, periodic and confined BCs
preemptively predict conduction block for high proportions of obstruction.
Linear BCs instead predict conduction even in the case where the finescale
simulation blocks.

### Volume averaging allows prediction of source/sink mismatch events

A critical component of the pro-arrhythmic effects of fibrosis is so-called
source/sink mismatch, in which spatial variation in the amount of excitable
tissue can create structures that permit conduction in one direction and not
another [2]. Unlike the other homogenised models that have been used to
represent the impacts of obstacles in cardiac electrophysiology [18, 19, 20],
the explicit representation of the local proportion of conductive tissue
($\phi$) in the homogenised models we construct provides them the potential to
capture this important electrophysiological dynamic. We examine this potential
using a set of small-scale nozzle-like structures (Fig 6a) that produce a
delay in activation, or outright block, when activation reaches the end of the
structure and attempts to emerge out into unobstructed tissue. Propagation
success is seen to depend predominantly on the width of the exit opening, $r$,
while the width of the entrance, $l$, has an effect only on borderline cases.
Both widths together control the extent of activation delay (calculated by
comparing the activation time at the opposite end of the tissue to the
activation time when no obstacle is present). We explore how well our
homogenised models predict these dynamics, and highlight the fact that the
sizes of the averaging volumes trialled are certainly large enough to obscure
the fine detail of the structure, and in most cases alter the effective width
of the entrance and exit.

a)

b)

c)

d)

Figure 6:  Homogenisation for the capture of source/sink mismatch.  a) The
nozzle structure used to evaluate the performance of homogenisation. b)
Activation delay and conduction block as predicted by homogenised models using
linear subproblem boundary conditions, for a range of exit widths $r$
($l=$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$$). Performance is good for
homogenised models with grid spacing $\Delta
x\leq$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. c) Performance of the
homogenised model using linear boundary conditions and $\Delta
x=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ for multiple combinations of
$l$ and $r$. Delay is well estimated by the homogenised model, and block is
perfectly predicted in all trialled scenarios. d) Example activation maps for
two finescale problems and the equivalent homogenised models (as in c)).
Patterns of hastened (blue) and delayed (red) activation are well recovered,
as is the occurrence of activation failure (light grey).

Homogenisation performance depends strongly on both the choice of boundary
conditions, and the size of the averaging volume. Homogenisation using
periodic or confined boundary conditions was found to perform poorly overall
(Figures S1-S2), most likely due to the poor handling of diagonal transport
seen for these choices in our test case (Fig 3). The nozzle structure
considered here also clearly violates the assumption of periodicity. Linear
boundary conditions, however, perform well providing the averaging volume is
not made too large (Fig 6b). In particular, homogenisation by a factor of ten
(resulting in a mesh spacing of $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
consistent with high-fidelity anatomic meshes) proved capable of perfectly
predicting the success or failure of propagation for all combinations of $l$
and $r$ values trialled (Fig 6c), and predicted delay accurately in the
majority of cases.

Two example maps of the activation delay using this best-performing
homogenisation are presented in Fig 6d. Here it is seen that although the
precise shape of the nozzle structure is well and truly lost in the
homogenised model, the pattern and timing of activation is well recovered.
Major discrepancies with the fine-scale problem occur only where a conductive
region is spuriously created by the homogenisation process (as any averaging
volume containing even a small amount of conductive material will become a
non-occupied element in the homogenised model). However, through its choice of
effective tensors and incorporation of the volume fraction, the homogenised
model is seen to compensate for this effect. In the case of a widening
channel, these locations lag in activation, resulting in additional sink to
slow propagation. This phenomenologically captures the increasing sink
experienced by a wavefront travelling through a widening structure. In the
opposite case of a narrowing channel, these regions activate more rapidly and
thus replicate the source-favoured balance for a wavefront moving through a
narrowing structure.

### Spiral wave anchoring can be predicted by homogenised models

A primary cause of arrhythmia are spiral waves, where cardiac tissue falls
into a self-sustaining pattern of continuous re-activation. Depending on both
electrophysiological and structural conditions, these spirals may stay fixed
rotating about a single region of the tissue or wander about it (with a chance
of self-annihilation upon collision with a tissue boundary). Fibrosis can act
to stabilise these spiral waves, with even small amounts of diffusely placed
obstacles shown to reduce spiral core wander [43]. Larger non-excitable
obstacles can also act as fixed locations to which a spiral wave anchors, and
is thus more likely to persist [44]. We explore here whether a homogenised
model can still produce this latter important component of fibrosis’ pro-
arrhythmic effect.

Spiral waves were simulated in two-dimensional, $6\text{\,}\mathrm{cm}$
$\\!\times\\!$ $6\text{\,}\mathrm{cm}$ slices of tissue, with a “finescale”
grid spacing of $\Delta x=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ used
for reasons of computational cost. Spiral waves were initiated using the
common cross field stimulus protocol [45], with a travelling wave initiated at
one edge of a two-dimensional slice of tissue, and then a second stimulus
triggered in one quadrant of the domain, timed to coincide with the
repolarisation front from the first wave. We use the “steep restitution” set
of parameter values provided by ten Tusscher et al. [25], which cause spiral
waves to break up and devolve into irregular patterns of activation. However,
then including a region of fibrotic obstruction causes the spiral wave to
anchor and be sustained apparently indefinitely without wave breakup (Fig 7,
first row).

Figure 7: Anchoring of spiral waves in homogenised models, as shown by
simulation snapshots. A region of 60% obstruction causes a spiral wave that
would otherwise break up to anchor and persist indefinitely. Homogenised
models successfully predict this anchoring despite losing the detailed
structure of the obstructed region. One homogenised model exhibits transient
breakthroughs from the obstructed region not seen in the finescale equivalent
(light blue circle). Videos are provided in the supplementary material.

Following their good performance on the previous scenarios considered, we
created homogenised models using linear boundary conditions, and averaging
volumes of size 10$\\!\times\\!$10 and 25$\\!\times\\!$25 elements. As the
base grid in these simulations is already
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$, we avoid the numerical effects of
significantly larger grid spacings by retaining the base grid but overlaying
the tensors obtained by block homogenisation. This approach allows for
feasible simulation of sufficient tissue to support spiral waves on the
finescale, to which the predictions of homogenised models can be directly
compared. We stress, however, that the most practical use case of
homogenisation remains incorporating very fine scale ($\sim\\!$
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$) structures into typical cardiac
meshes ($\sim\\!$ 100-$250\text{\,}\mathrm{\SIUnitSymbolMicro m}$).

Most immediately, the homogenised models are seen to correctly predict the
anchoring of spiral waves to a region of obstruction (Fig 7). Although precise
positions of wavefronts do not match in individual snapshots, frequencies of
rotation measured over multiple rotations are very similar
($3.2\text{\,}\mathrm{Hz}$ for $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$
finescale model and $1000\text{\,}\mathrm{\SIUnitSymbolMicro m}$ homogenised
models and $3.1\text{\,}\mathrm{Hz}$ for
$2500\text{\,}\mathrm{\SIUnitSymbolMicro m}$ homogenised model as estimated
from resultant movies, available in the supplement). A noticeable discrepancy
between finescale and homogenised models can be seen for the
$1000\text{\,}\mathrm{\SIUnitSymbolMicro m}$ model, where the obstructed
region proves too conductive and hence admits breakthroughs that are not seen
in the finescale model due to the difference in timing caused by slower
traversal through the obstructed region. However, this does not compromise the
core behaviours of spiral wave anchoring and frequency of rotation that we
want the homogenised models to predict.

## Conclusions

Homogenisation has seen only limited use in the modelling of fibrosis in
cardiac electrophysiology, despite the technique presenting a natural means
for incorporating the effects of sub-mesh-scale obstructions such as cardiac
fibrosis into existing heart meshes. The dynamic behaviours seen in
electrophysiological models, owing to their highly influential and strongly
non-linear reaction terms, do present a significant challenge to prototypical
homogenisation, which concerns itself only with the calculation of macroscopic
transport properties. However, by careful application of the volume averaging
theory for non-conductive obstructions [27], and a thorough consideration of
how different choices of closure problem boundary conditions behave in several
pernicious cases of interest in this field, we have demonstrated how
homogenisation can robustly capture some of the key dynamics promoting
arrhythmia.

To demonstrate the method, we have used regular, fine-scale
($10\text{\,}\mathrm{\SIUnitSymbolMicro m}$) meshes that captured the spatial
scale of fibrotic obstructions to conduction, in two dimensions. In practice,
a modeller will most likely have a coarser-scale
($100-$250\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) mesh, regular or irregular
and in two or three dimensions, on which they wish to simulate cardiac
activity. In this case, our suggested homogenisation approach would proceed by
first temporarily creating a finer-scale grid on which to represent fibrosis,
and then solving the set of closure subproblems Eq (10) on this scale (each
consisting of a comparatively small number of nodes and elements and able to
be processed in parallel). The derived effective conductivities can then be
used in the original coarser mesh, along with the volume fractions, through Eq
(7). Although some of the boundary conditions we have considered here
(confined and periodic) do not naturally generalise to irregular grids, the
fine-scale grid used to represent fibrosis can be freely chosen. As such, a
compelling option is to use a regular grid on the finescale, such that
boundary conditions and the concept of periodicity can be treated in the same
fashion we have demonstrated here. The portions of the finescale mesh in which
the irregularly-shaped elements of the coarser mesh are embedded then becomes
the “skin” that we have demonstrated here to reduce boundary effects.

Overall, homogenisation has been shown here to perform better than might be
expected, considering the reaction-dominated dynamics of excitation
propagation. Key pro-arrhythmic impacts of fibrotic obstacles to conduction
(such as collagenous deposits) have been observed in appropriately homogenised
models. This included slowed conduction in fibrosis-afflicted tissue,
conduction delay and/or block due to events of source/sink mismatch, and
anchoring of spiral waves to heavily obstructed regions. In particular,
predictions of conduction slowing were very good for levels of obstruction up
to 30-40%, matching the proportions of collagen that can be seen in typical
histological sections of fibrotic tissue [6].

The best-performing homogenised models here used reductions in node count by
factors of 100 and 625 compared to simulations explicitly resolving finescale
detail, representing significant computational speedup. In three dimensions,
this speedup is expected to be even more pronounced. Interpreted differently,
the homogenisation approach detailed here makes it feasible to incorporate
finescale structures into existing heart meshes already composed of large
numbers of nodes and elements, where refinement of the mesh down to the
finescale is infeasible due to both computational speed and memory
limitations. Although there is a computational cost associated with solving
the closure problems that define effective conductivities in a homogenised
model, this is a one-off cost, and for the problems considered here this cost
was only minor relative to the cost of simulating on the finescale. Where the
solution of closure problems threatens to become a bottleneck, they may be
solved in parallel or via semi-analytical techniques that further reduce the
time required [46].

We have considered the performance of different choices for closure subproblem
boundary conditions, both inclusive of, and separate from, the well-known
issue of the sensitivity of the monodomain (and bidomain) model to the spatial
discretisation [5]. Despite their reputation for reduced accuracy [33], linear
boundary conditions proved most accurate and robust overall, here assisted by
the inclusion of skin in closure subproblems [35] that improved their
performance (Fig 3). Even when periodic boundary conditions are applicable and
perform best in terms of homogenisation error, linear boundary conditions
become competitive again when the effects of changing gridsize are also
incorporated (Fig 4). The key weakness of linear boundary conditions in this
context is the potential for complete barriers to conduction to “leak” (Fig
2), although we consider this superior to the other choices for boundary
conditions that have an even higher tendency to over-predict conduction block.
Linear boundary conditions were also seen to be clearly superior for the
nozzle problem we used to explore source/sink mismatch.

The primary limitation of homogenisation in this context is of course the
elimination of the microscopic structure by the homogenisation process.
Although we have demonstrated capture of source/sink mismatch and spiral wave
anchorage through proof-of-concept experiments, we have not exhaustively
tested whether homogenised models can predict the precise manifestations of
these effects across the many different types of obstacle arrangement that may
be of interest. The further challenge of capturing re-entries that “live” on
the micro-scale [3] has also not been considered here.

In summary, we have demonstrated that when used with care, homogenisation by
volume averaging has good potential even for the case where the sensitive and
reaction-dominated dynamics of cardiac electrophysiology meet wholly non-
conductive material in non-periodic arrangements. Following the result that
discrepancies between fine-scale and homogenised models owed mostly to the
numerical consequences of the different length scales of the models, this
suggests that such homogenisation would perform even better in other settings
where the governing equations produce travelling waves with less sharp fronts.
The core dynamics of excitation propagation and refractoriness can also be
seen, for example, in Ca2+ [47] or cyclic adenosine monophosphate [48]
signalling in biological systems. The most critical targets for future work
are numerical schemes for the monodomain and bidomain equations with further
reduced grid sensitivity [5, 49], and/or a new homogenisation paradigm that
somehow respects the distinct, tortuous paths through highly obstructed tissue
and hence might be capable of generating micro re-entries even in larger-scale
homogenised models.

## Acknowledgments

BAJL, IT, PB and KB are funded by the Australian Research Council Centre of
Excellence for Mathematical and Statistical Frontiers (project number
CE140100049), an initiative of the Australian government. RWS acknowledges the
support by the Brazilian Government via CAPES, CNPq, FAPEMIG, and Universidade
Federal de Juiz de Fora, and by the Australian Government via the Endeavour
Research Leadership Award from the Department of Education. ABO acknowledges a
British Heart Foundation Intermediate Basic Science Fellowship
(FS/17/22/32644), and an Impact for Infrastructure Award from the National
Centre for the Replacement, Refinement and Reduction of Animals in Research
(NC/P001076/1).

## Code Availability

MATLAB code for the simulation of the monodomain model and the formulation of
homogenised models through homogenisation by volume averaging will be made
available via GitHub at https://github.com/betalawson/fibro-homogenisation.

## References

  * [1] X. Zhou, A. Bueno-Orovio, and B. Rodriguez, “In silico evaluation of arrhythmia,” Curr. Opin. Physiol., vol. 1, pp. 95–103, 2018.
  * [2] T. P. Nguyen, Z. Qu, and J. N. Weiss, “Cardiac fibrosis and arrhythmogenesis: the road to repair is paved with perils,” J. Mol. Cell Cardiol., vol. 70, pp. 83–91, 2014.
  * [3] B. J. Hansen, J. Zhao, T. A. Csepe, B. T. Moore, N. Li, L. A. Jayne, A. Kalyanasundaram, P. Lim, A. Bratasz, K. A. Powell, O. P. Simonetti, R. S. D. Higgins, A. Kilic, P. J. Mohler, P. M. L. Janssen, R. Weiss, J. D. Hummel, and V. V. Fedorov, “Atrial fibrillation driven by micro-anatomic intramural re-entry revealed by simultaneous sub-epicardial and sub-endocardial optical mapping in explanted human hearts,” Eur. Heart J., vol. 36, pp. 2390–2401, 2015.
  * [4] S. Zahid, H. Cochet, P. M. Boyle, E. L. Schwarz, K. N. Whyte, E. J. Vigmond, R. Duboi, M. Hocini, M. Haïssaguerre, P. Jaïs, and N. A. Trayanova, “Patient-derived models link re-entrant driver localization in atrial fibrillation to fibrosis spatial pattern,” Cardiovasc. Res., vol. 110, pp. 443–454, 2016.
  * [5] P. Pathmanathan, M. O. Bernabeu, S. A. Niederer, D. J. Gavaghan, and D. Kay, “Computational modelling of cardiac electrophysiology: explanation of the variability of results from different numerical solvers,” Int. J. Numer. Method Biomed. Eng., vol. 28, pp. 890–903, 2012.
  * [6] S. de Jong, T. A. B. van Veen, H. V. M. van Rijen, and J. M. T. de Bakker, “Fibrosis and cardiac arrhythmias,” J. Cardiovasc. Pharmacol., vol. 57, pp. 630–638, 2011.
  * [7] J. P. Whiteley, “An efficient numerical technique for the solution of the monodomain and bidomain equations,” IEEE Trans. Biomed. Eng., vol. 53, pp. 2139–2147, 2006.
  * [8] R. Sachetto Oliveira, M. Rocha, Bernardo, D. Burgarelli, W. Meira Jr., C. Constantinides, and R. Weber dos Santos, “Performance evaluation of GPU parallelization, space-time adaptive algorithms, and their combination for simulating cardiac electrophysiology,” Int. J. Numer. Meth. Biomed. Engng., vol. 34, p. c2913, 2018.
  * [9] C. M. Costa, F. O. Campos, A. J. Prassl, R. W. dos Santos, D. Sánchez-Quintana, H. Ahammer, E. Hofer, and G. Plank, “An efficient finite element approach for modeling fibrotic clefts in the heart,” IEEE Trans. Biomed. Eng., vol. 61, pp. 900–910, 2014.
  * [10] Y. Davit, C. G. Bell, H. M. Byrne, L. A. C. Chapman, L. S. Kimpton, G. E. Lang, K. H. L. Leonard, J. M. Oliver, N. C. Pearson, R. J. Shipley, S. L. Waters, J. P. Whiteley, B. D. Wood, and M. Quintard, “Homogenization via formal multiscale asymptotics and volume averaging: How do the two techniques compare?,” Adv. Water Res., vol. 62B, pp. 178–206, 2013.
  * [11] J. Keener and J. Sneyd, Mathematical Physiology I: Cellular Physiology.  Springer-Verlag New York, 2009.
  * [12] P. E. Hand, B. E. Griffith, and C. S. Peskin, “Deriving macroscopic myocardial conductivities by homogenization of microscopic models,” Bull. Math. Biol., vol. 71, pp. 1707–1726, 2009.
  * [13] C. M. Costa and dos Santos Rodrigo W., “Limitations of the homogenized cardiac Monodomain model for the case of low gap junctional coupling,” Conf. Proc. IEEE Eng. Med. Biol. Soc., vol. 2010, pp. 228–231, 2010.
  * [14] G. Richardson and S. J. Chapman, “Derivation of the bidomain equations for a beating heart with a general microstructure,” SIAM J. Appl. Math., vol. 71, pp. 657–675, 2011.
  * [15] E. Grandelius and K. H. Karlsen, “The cardiac bidomain model and homogenization,” Netw. Heterog. Media, vol. 14, pp. 173–204, 2019.
  * [16] J. M. Kim, N. Bursac, and C. S. Henriquez, “A computer model of engineered cardiac monolayers,” Biophys. J., vol. 98, pp. 1762–1771, 2010.
  * [17] D. E. Hurtado, J. Jilberto, and G. Panasenko, “Non-ohmic tissue conduction in cardiac electrophysiology: Upscaling the non-linear voltage-dependent conductance of gap junctions,” PLoS Comp. Biol., vol. 16, p. e1007232, 2020\.
  * [18] A. Davidović, Y. Coudière, and Y. Bourgault, “Image-based modeling of the heterogeneity of propagation of the cardiac action potential. example of rat heart high resolution mri,” in Functional Imaging and Modelling of the Heart (M. Pop and G. A. Wright, eds.), (Cham), pp. 260–270, Springer International Publishing, 2017.
  * [19] A. Gokhale, Tanmay, H. Asfour, S. Verma, N. Bursac, and C. S. Henriquez, “Microheterogeneity-induced conduction slowing and wavefront collisions govern macroscopic conduction behaviour: A computational and experimental study,” PLoS Comp. Biol., vol. 14, p. e1006276, 2018.
  * [20] T. Austin, M. Trew, and A. Pullan, “Multilevel homogenization applied to the cardiac bidomain equations,” in Proceedings of the 28th IEEE EMBS Annual International Conference, pp. 584–587, 2006.
  * [21] J. Sundnes, G. T. Lines, X. Cai, B. F. Nielsen, K. Mardal, and A. Tveito, Computing the Electrical Activity in the Heart.  Springer-Verlag, 2006.
  * [22] M. Potse, B. Dubé, J. Richer, A. Vinet, and R. M. Gulrajani, “A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart,” IEEE Trans. Biomed. Eng., vol. 53, pp. 2425–2435, 2006.
  * [23] J. Sundnes, B. F. Nielsen, K. A. Mardal, X. Cai, G. T. Lines, and A. Tveito, “On the computational complexity of the bidomain and the monodomain models of electrophysiology,” Ann. Biomed. Eng., vol. 34, pp. 1088–1097, 2006\.
  * [24] Y. Bourgault and C. Pierre, “Comparing the bidomain and monodomain models in electro-cardiology through convergence analysis,” HAL, vol. 00545888, 2010\.
  * [25] K. H. W. J. Ten Tusscher and A. V. Panfilov, “Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions,” Phys. Med. Biol., vol. 51, pp. 6141–6156, 2006\.
  * [26] M. Perego and A. Veneziani, “An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations,” Elecron. Trans. Numer. Anal., vol. 35, pp. 234–256, 2009.
  * [27] S. Whitaker, The Method of Volume Averaging.  Springer Science+Business Media, 1999.
  * [28] F. H. Fenton, M. Cherry, Elizabeth, A. Karma, and W.-J. Rappel, “Modeling wave propagation in realistic heart geometries using the phase-field method,” Chaos, vol. 15, p. 013502, 2005.
  * [29] A. Bueno-Orovio, V. M. Pérez-García, and F. H. Fenton, “Spectral methods for partial differential equations in irregular domains: The spectral smoothed boundary method,” SIAM J. Sci. Comput., vol. 28, pp. 886–900, 2006.
  * [30] X. H. Wu, Y. Efendiev, and T. Y. Hou, “Analysis of upscaling absolute permeability,” Discrete Continuous Dyn. Syst. Ser. B, vol. 2, pp. 185–204, 2002.
  * [31] L. J. Durlofsky, “Representation of grid block permeability in coarse scale models of randomly heterogeneous porous media,” Water Resour. Res., vol. 28, pp. 1791–1800, 1992.
  * [32] S. J. Asirvatham and W. G. Stevenson, “Editor’s perspective: The isthmus of uncertainty,” Circ. Arrhythm. Electrophysiol., vol. 7, pp. 175–177, 2014\.
  * [33] A. Szymkiewicz, Modelling Water Flow in Unsaturated Porous Media: Accounting for Nonlinear Permeability and Material Heterogeneity.  Springer-Verlag Berlin Heidelberg, 2013.
  * [34] X. Sánchez-Vila, J. P. Girardi, and J. carrera, “A synthesis of approaches to upscaling of hydraulic conductivities,” Water Resour. Res., vol. 31, pp. 867–882, 1995.
  * [35] J. J. Gómez-Hernández, “Simulation of block permeabilities conditioned upon data measured at a different scale,” in ModelCARE 90: Calibration and Reliability in Groundwater Modelling (K. Kovar, ed.), pp. 407–416, IAHS, 1990.
  * [36] X.-H. Wen, L. J. Durlofsky, S. H. Lee, and M. G. Edwards, “Full tensor upscaling of geologically complex resevoir descriptions,” Conf. Proc. Soc. Pet. Eng., vol. 62928, 2000.
  * [37] P. Henning and D. Peterseim, “Oversampling for the multiscale finite element method,” Multiscale Model. Simul., vol. 11, pp. 1149–1175, 2013.
  * [38] N. J. Higham, “Computing a nearest symmetric positive semidefinite matrix,” Linear Algebra Its Appl., vol. 103, pp. 103–118, 1988.
  * [39] P. Renard and G. de Marsily, “Calculating equivalent permeability: a review,” Adv. Water Res., vol. 20, pp. 253–278, 1997.
  * [40] J. M. T. de Bakker, F. J. L. van Capelle, M. J. Janse, S. Tasseron, J. T. Vermeulen, N. de Jonge, and J. R. Lahpor, “Slow conduction in the infarcted human heart: ‘zigzag’ course of activation,” Circulation, vol. 88, pp. 915–926, 1993.
  * [41] J. L. R. M. Smeets, M. A. Allessie, W. J. E. P. Lammers, F. I. M. Bonke, and J. Hollen, “The wavelength of the cardiac impulse and reentrant arrhythmias in isolated rabbit atrium. the role of heart rate, autonomic transmitters, temperature and potassium,” Circ. Res., vol. 58, pp. 96–108, 1986.
  * [42] T. Nezlobinsky, O. Solovyova, and A. V. Panfilov, “Anisotropic conduction in the myocardium due to fibrosis: the effect of texture on wave propagation,” Sci. Rep., vol. 10, p. 764, 2020.
  * [43] K. H. W. J. Ten Tusscher and A. V. Panfilov, “Influence of diffuse fibrosis on wave propagation in human ventricular tissue,” Europace, vol. 9, pp. vi38–45, 2007.
  * [44] J. M. Davidenko, A. V. Pertsov, R. Salomonsz, W. Baxter, and J. Jalife, “Stationary and drifting spiral waves of excitation in isolated cardiac muscle,” Nature, vol. 355, pp. 349–351, 1992.
  * [45] J. Beaumont, N. Davidenko, J. M. Davidenko, and J. Jalife, “Spiral waves in two-dimensional models of ventricular muscle: Formation of a stationary core,” Biophys. J., vol. 75, pp. 1–14, 1998.
  * [46] N. G. March, E. J. Carr, and I. W. Turner, “Fast semi-analytical solution of the homogenization boundary value problem for block locally-isotropic heterogeneous media,” Appl. Math. Model., p. https://doi.org/10.1016/j.apm.2020.09.022, 2020.
  * [47] J. Lechleiter, S. Girard, E. Peralta, and D. Clapham, “Sprial calcium wave propagation and annihilation in xenopus laevis oocytes,” Science, vol. 252, pp. 123–126, 1991.
  * [48] D. Dormann, J.-Y. Kim, P. N. Devreotes, and C. J. Weijer, “cAMP receptor affinity controls wave dynamics, geometry and morphogenesis in dictyostelium,” J. Cell Sci., vol. 114, pp. 2513–2523, 2001.
  * [49] C. M. Costa, P. A. A. Silva, and R. W. dos Santos, “Mind the gap: A semicontinuum model for discrete electrical propagation in cardiac tissue,” IEEE Trans. Biomed. Eng., vol. 63, pp. 765–774, 2016.

## Supporting information

a)

b)

Figure S1:  The performance of confined boundary conditions for capture of the
effects of source/sink mismatch. a) Activation delay and conduction block as
predicted by homogenised models of different characteristic lengths, for a
range of exit widths $r$ ($l=$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$$).
Large characteristic lengths for homogenised models result in significant
overestimation of delay and block, while the smaller characteristic lengths
result in underestimation. b) Activation delay (AD) and conduction block as
predicted by the best-performing homogenised model ($\Delta
x=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$). The rapid transition from low
AD to block ($r\sim$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ in the
finescale model) is predicted poorly by even this best-performing homogenised
model, indicating a failure of confined boundary conditions.

a)

b)

Figure S2:  The performance of periodic boundary conditions for capture of the
effects of source/sink mismatch. a) Activation delay and conduction block as
predicted by homogenised models of different characteristic lengths, for a
range of exit widths $r$ ($l=$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$$).
Large characteristic lengths for homogenised models result in significant
overestimation of delay and block, while the smaller characteristic lengths
result in underestimation. b) Activation delay (AD) and conduction block as
predicted by the best-performing homogenised model ($\Delta
x=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$). The rapid transition from low
AD to block ($r\sim$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ in the
finescale model) is predicted poorly by even this best-performing homogenised
model, indicating a failure of periodic boundary conditions.


* [87] P. A. Hill, P. A. Nelson, O. Kirkeby, and H. Hamada, “Resolution of Front–Back Confusion in Virtual Acoustic Imaging Systems,” _Journal of the Acoustical Society of America_ , vol. 108, no. 6, pp. 2901–2910, 2000.
  * [88] U.-H. Kim, K. Nakadai, and H. G. Okuno, “Improved Sound Source Localization and Front-Back Disambiguation for Humanoid Robots with Two Ears,” in _Recent Trends in Applied Artificial Intelligence_ , M. Ali, T. Bosse, K. V. Hindriks, M. Hoogendoorn, C. M. Jonker, and J. Treur, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2013, pp. 282–291.
  * [89] W. Bangs and P. Schultheis, “Space-Time Processing for Optimal Parameter Estimation,” _Signal Processing_ , pp. 577–590, 1973.
  * [90] W. Hahn and S. Tretter, “Optimum Processing for Delay-Vector Estimation in Passive Signal Arrays,” _IEEE Trans. on Information Theory_ , vol. 19, no. 5, pp. 608–614, Sep. 1973.
  * [91] M. Wax and T. Kailath, “Optimum Localization of Multiple Sources by Passive Arrays,” _IEEE Trans. on Acoustics, Speech, and Signal Processing_ , vol. 31, no. 5, pp. 1210–1217, Oct. 1983.
  * [92] M. Taseska and E. A. P. Habets, “Spotforming: Spatial Filtering With Distributed Arrays for Position-Selective Sound Acquisition,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 24, no. 7, pp. 1291–1304, Jul. 2016.
  * [93] M. Omologo, P. G. Svaizer, and R. D. Mori, _Acoustic Transduction_. San Diego, CA: Academic Press, 1998, pp. 23–69.
  * [94] H. F. Silverman, Y. Yu, J. M. Sachar, and W. R. Patterson, “Performance of Real-Time Source-Location Estimators for a Large-Aperture Microphone Array,” _IEEE Trans. on Acoustics, Speech, and Signal Processsing_, vol. 13, no. 4, pp. 593–606, Jul. 2005.
  * [95] D. Salvati, C. Drioli, and G. L. Foresti, “A Low-Complexity Robust Beamforming Using Diagonal Unloading for Acoustic Source Localization,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 26, no. 3, pp. 609–622, Mar. 2018.
  * [96] B. Rafaely, _Fundamentals of Spherical Array Processing_ , ser. Springer Topics in Signal Processing. Springer, 2015\.
  * [97] ——, “Analysis and Design of Spherical Microphone Arrays,” _IEEE Trans. on Speech and Audio Processing_ , vol. 13, no. 1, pp. 135–143, Jan. 2005\.
  * [98] ——, “Plane-wave Decomposition of the Sound Field on a Sphere by Spherical Convolution,” _Journal of the Acoustical Society of America_ , vol. 116, no. 4, pp. 2149–2157, 2004.
  * [99] L. Kumar and R. M. Hegde, “Near-Field Acoustic Source Localization and Beamforming in Spherical Harmonics Domain,” _IEEE Trans. on Signal Processing_, vol. 64, no. 13, pp. 3351–3361, Jul. 2016.
  * [100] D. Jarrett, E. A. P. Habets, and P. A. Naylor, _Theory and Applications of Spherical Microphone Array Processing_. Springer, 2016.
  * [101] A. H. Moore, C. Evers, and P. A. Naylor, “Direction of Arrival Estimation in the Spherical Harmonic Domain Using Subspace Pseudointensity Vectors,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 25, no. 1, pp. 178–192, Jan. 2017.
  * [102] H. Sawada, R. Mukai, and S. Makino, “Direction of Arrival Estimation for Multiple Source Signals using Independent Component Analysis,” in _Proc. of International Symposium on Signal Processing and Its Applications_ , vol. 2, Jul. 2003, pp. 411–414 vol.2.
  * [103] A. Lombard, Y. Zheng, H. Buchner, and W. Kellermann, “TDOA Estimation for Multiple Sound Sources in Noisy and Reverberant Environments Using Broadband Independent Component Analysis,” _IEEE Trans. on Audio, Speech, and Language Processing_, vol. 19, no. 6, pp. 1490–1503, Aug. 2011.
  * [104] M. I. Mandel, R. J. Weiss, and D. P. W. Ellis, “Model-Based Expectation-Maximization Source Separation and Localization,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 18, no. 2, pp. 382–394, Feb. 2010.
  * [105] R. J. Weiss, M. I. Mandel, and D. P. Ellis, “Combining Localization Cues and Source Model Constraints for Binaural Source Separation,” _Speech Communication_ , vol. 53, no. 5, pp. 606 – 621, 2011.
  * [106] R. Schmidt, “Multiple Emitter Location and Signal Parameter Estimation,” _IEEE Trans. on Antennas and Propagation_ , vol. 34, no. 3, pp. 276–280, Mar. 1986.
  * [107] R. Roy and T. Kailath, “ESPRIT-Estimation of Signal Parameters via Rotational Invariance Techniques,” _IEEE Trans. on Acoustics, Speech, and Signal Processing_ , vol. 37, no. 7, pp. 984–995, Jul. 1989.
  * [108] H. Teutsch and W. Kellermann, “EB-ESPRIT: 2D Localization of Multiple Wideband Acoustic Sources using Eigen-beams,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_, Mar. 2005.
  * [109] ——, “Detection and Localization of Multiple Wideband Acoustic Sources Based on Wavefield Decomposition using Spherical Apertures,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_, Las Vegas (Nevada), USA, Mar. 2008, pp. 5276–5279.
  * [110] G. Golub and W. Kahan, “Calculating the Singular Values and Pseudo-Inverse of a Matrix,” _Journal of the Society for Industrial and Applied Mathematics Series B Numerical Analysis_ , vol. 2, no. 2, pp. 205–224, 1965.
  * [111] H. L. Van Trees, _Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory_. New York: Wiley, 2002.
  * [112] O. Nadiri and B. Rafaely, “Localization of Multiple Speakers under High Reverberation Using a Spherical Microphone Array and the Direct-Path Dominance Test,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 22, no. 10, Oct. 2014.
  * [113] B. Ristic, _Particle Filters for Random Set Models_ , 1st ed. New York: Springer, 2013.
  * [114] R. P. S. Mahler, _Statistical Multisource Multitarget Information Fusion_. Artech House, 2007.
  * [115] A. Deleforge and R. Horaud, “2D Sound-Source Localization on the Binaural Manifold,” in _IEEE Workshop on Machine Learning for Signal Processing (MLSP)_. Santander, Spain: IEEE, Sep. 2012, pp. 1–6.
  * [116] B. Laufer-Goldshtein, R. Talmon, and S. Gannot, “Semi-Supervised Sound Source Localization Based on Manifold Regularization,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 24, no. 8, pp. 1393–1407, Aug. 2016.
  * [117] X. Xiao, S. Zhao, X. Zhong, D. L. Jones, E. S. Chng, and H. Li, “A Learning-based Approach to Direction of Arrival Estimation in Noisy and Reverberant Environments,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_, Apr. 2015, pp. 2814–2818.
  * [118] E. L. Ferguson, S. B. Williams, and C. T. Jin, “Sound Source Localization in a Multipath Environment Using Convolutional Neural Networks,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_, Apr. 2018, pp. 2386–2390.
  * [119] R. Takeda and K. Komatani, “Sound Source Localization based on Deep Neural Networks with Directional Activate Function Exploiting Phase Information,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_, Mar. 2016, pp. 405–409.
  * [120] ——, “Unsupervised Adaptation of Deep Neural Networks for Sound Source Localization using Entropy Minimization,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_, Mar. 2017, pp. 2217–2221.
  * [121] N. Ma, T. May, and G. J. Brown, “Exploiting Deep Neural Networks and Head Movements for Robust Binaural Localization of Multiple Sources in Reverberant Environments,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 25, no. 12, pp. 2444–2453, Dec. 2017.
  * [122] M. R. Bai, S. Lan, and J. Huang, “Time Difference of Arrival (TDOA)-Based Acoustic Source Localization and Signal Extraction for Intelligent Audio Classification,” in _Proc. IEEE Sensor Array and Multichannel Signal Processing Workshop (SAM)_ , Jul. 2018, pp. 632–636.
  * [123] F. Grondin, F. Glass, I. Sobieraj, and M. D. Plumbley, “Sound Event Localization and Detection Using CRNN on Pairs of Microphones,” in _The Detection and Classification of Acoustic Scenes and Events Workshop (DCASE)_ , New York, USA, Oct. 2019.
  * [124] N. Ma, J. A. Gonzalez, and G. J. Brown, “Robust Binaural Localization of a Target Sound Source by Combining Spectral Source Models and Deep Neural Networks,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 26, no. 11, pp. 2122–2131, Nov. 2018.
  * [125] S. Chakrabarty and E. A. P. Habets, “Multi-Speaker Localization using Convolutional Neural Network Trained with Noise,” in _Proc. of ML4Audio Workshop at NIPS_ , 2017.
  * [126] D. Hinrichsen and A. J. Pritchard, _Mathematical Systems Theory I: Modelling, State Space Analysis, Stability and Robustness_ , ser. Texts in Applied Mathematics. Berlin Heidelberg: Springer, 2005.
  * [127] R. E. Kalman, “A New Approach to Linear Filtering and Prediction Problems,” _Transactions of the ASME–Journal of Basic Engineering_ , vol. 82, no. Series D, pp. 35–45, 1960.
  * [128] B. Ristic, S. Arulampalam, and N. Gordon, _Beyond the Kalman filter: Particle Filters for Tracking Applications_. Boston: Artech House, 2004.
  * [129] D. B. Ward, E. A. Lehmann, and R. C. Williamson, “Particle Filtering Algorithms for Tracking an Acoustic Source in a Reverberant Environment,” _IEEE Trans. on Speech and Audio Processing_, vol. 11, no. 6, pp. 826–836, Nov. 2003.
  * [130] E. A. Lehmann and R. C. Williamson, “Particle Filter Design Using Importance Sampling for Acoustic Source Localisation and Tracking in Reverberant Environments,” _EURASIP Journal on Advances in Signal Processing_ , Jun. 2006.
  * [131] A. Doucet, N. de Freitas, and N. Gordon, _Sequential Monte Carlo Methods in Practice_ , ser. Information Science and Statistics. Springer, 2001.
  * [132] A. Doucet, S. Godsill, and C. Andrieu, “On Sequential Monte Carlo Sampling Methods for Bayesian Filtering,” _Statistics and Computing_ , vol. 10, no. 3, 2000.
  * [133] M. F. Fallon and S. J. Godsill, “Acoustic Source Localization and Tracking of a Time-Varying Number of Speakers,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 20, no. 4, pp. 1409–1415, May 2012.
  * [134] P. Li, R. Goodall, and V. Kadirkamanathan, “Estimation of Parameters in a Linear State Space Model using a Rao-Blackwellised Particle Filter,” _IEE Proceedings - Control Theory and Applications_ , vol. 151, no. 6, pp. 727–738, Nov. 2004.
  * [135] X. Zhong and J. R. Hopgood, “Particle Filtering for TDOA based Acoustic Source Ttracking: Nonconcurrent Multiple Talkers,” _Signal Processing_ , vol. 96, pp. 382 – 394, 2014.
  * [136] A. Levy, S. Gannot, and E. A. P. Habets, “Multiple-Hypothesis Extended Particle Filter for Acoustic Source Localization in Reverberant Environments,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 19, no. 6, pp. 1540–1555, Aug. 2011.
  * [137] T. Li, M. Bolic, and P. M. Djuric, “Resampling Methods for Particle Filtering: Classification, Implementation, and Strategies,” _IEEE Signal Processing Magazine_ , vol. 32, no. 3, pp. 70–86, May 2015.
  * [138] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A Tutorial on Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking,” _IEEE Trans. on Signal Processing_, vol. 50, no. 2, pp. 174–188, Feb. 2002.
  * [139] M. Bolic and P. M. D. and, “New Resampling Algorithms for Particle Filters,” in _Proc. of IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_, vol. 2, Apr. 2003, pp. II–589.
  * [140] M. Halimeh, C. Hümmer, A. Brendel, and W. Kellermann, “Hybrid Particle Filtering based on an Elitist Resampling Scheme,” in _Sensor Array and Multichannel Signal Processing Workshop (SAM)_ , Jul. 2018, pp. 257–261.
  * [141] M. Halimeh, A. Brendel, and W. Kellermann, “Evolutionary Resampling for Multi-Target Tracking using Probability Hypothesis Density Filter,” in _Proc. of European Signal Processing Conf.(EUSIPCO)_, Sep. 2018, pp. 647–651.
  * [142] A. Plinge, F. Jacob, R. Haeb-Umbach, and G. A. Fink, “Acoustic Microphone Geometry Calibration: An Overview and Experimental Evaluation of State-of-the-Art Algorithms,” _IEEE Signal Processing Magazine_ , vol. 33, no. 4, pp. 14–29, Jul. 2016.
  * [143] D. Cherkassky and S. Gannot, “Blind Synchronization in Wireless Acoustic Sensor Networks,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 25, no. 3, pp. 651–661, Mar. 2017.
  * [144] D. J. Salmond, “Mixture Reduction Algorithms for Point and Extended Object Tracking in Clutter,” _IEEE Trans. on Aerospace and Electronic Systems_ , vol. 45, no. 2, pp. 667–686, Apr. 2009.
  * [145] K. V. Mardia and P. E. Jupp, _Directional Statistics_. John Wiley & Sons, 2009, vol. 494.
  * [146] K. V. Mardia, “Bayesian Analysis for Bivariate von Mises Distributions,” _Journal of Applied Statistics_ , vol. 37, no. 3, pp. 515–528, 2010.
  * [147] M. Jeub, C. Nelke, C. Beaugeant, and P. Vary, “Blind Estimation of the Coherent-to-Diffuse Energy Ratio from Noisy Speech Signals,” in _Proc. of European Signal Processing Conf.(EUSIPCO)_, Aug. 2011, pp. 1347–1351.
  * [148] S. Braun, A. Kuklasiński, O. Schwartz, O. Thiergart, E. A. P. Habets, S. Gannot, S. Doclo, and J. Jensen, “Evaluation and Comparison of Late Reverberation Power Spectral Density Estimators,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 26, no. 6, pp. 1056–1071, Jun. 2018.
  * [149] S. S. Blackman, “Association and Fusion of Multiple Sensor Data,” in _Multitarget-Multisensor Tracking: Advanced Algorithms_ , Y. Bar-Shalom, Ed. Norwood, MA: Artech House, 1990, ch. 6, pp. 187–218.
  * [150] E. C. Cherry, “Some Experiments on the Recognition of Speech, with One and with Two Ears,” _Journal of the Acoustical Society of America_ , vol. 25, no. 5, pp. 975–979, 1953.
  * [151] S. Haykin and Z. Chen, “The Cocktail Party Problem,” _Neural Computation_ , vol. 17, no. 9, pp. 1875–1902, 2005.
  * [152] Y. Bar-Shalom and E. Tse, “Tracking in a Cluttered Environment with Probabilistic Data Association,” _Automatica_ , vol. 11, no. 5, pp. 451 – 460, 1975.
  * [153] Y. Bar-Shalom, F. Daum, and J. Huang, “The Probabilistic Data Association Filter,” _IEEE Control Systems Magazine_ , vol. 29, no. 6, pp. 82–100, Dec. 2009.
  * [154] T. Fortmann, Y. Bar-Shalom, and M. Scheffe, “Sonar Tracking of Multiple Targets using Joint Probabilistic Data Association,” _IEEE Journal of Oceanic Engineering_ , vol. 8, no. 3, pp. 173–184, Jul. 1983.
  * [155] T. Gehrig and J. McDonough, “Tracking Multiple Speakers with Probabilistic Data Association Filters,” in _Multimodal Technologies for Perception of Humans_ , R. Stiefelhagen and J. Garofolo, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, pp. 137–150.
  * [156] Y. Ban, X. Alameda-Pineda, C. Evers, and R. Horaud, “Tracking multiple audio sources with the von mises distribution and variational em,” _IEEE Signal Processing Letters_ , vol. 26, no. 6, pp. 798–802, Jun. 2019\.
  * [157] R. P. S. Mahler, _Advances in Statistical Multisource-Multitarget Information Fusion_. Artech House, 2014\.
  * [158] ——, “Statistics 101 for Multisensor, Multitarget Data Fusion,” _Aerospace and Electronic Systems Magazine, IEEE_ , vol. 19, no. 1, Jan. 2004\.
  * [159] ——, ““Statistics 102” for Multisource-Multitarget Detection and Tracking,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 7, no. 3, pp. 376–389, Jun. 2013.
  * [160] W.-K. Ma, B.-N. Vo, S. S. Singh, and A. Baddeley, “Tracking an Unknown Time-Varying Number of Speakers Using TDOA Measurements: A Random Finite Set Approach,” _IEEE Trans. on Signal Processing_, vol. 54, no. 9, pp. 3291–3304, Sep. 2006.
  * [161] H. Pessentheiner, “Localization, Characterization, and Tracking of Harmonic Sources With Applications to Speech Signal Processing,” Ph.D. dissertation, Graz University of Technology, Jan. 2017.
  * [162] I. Marković, J. Ćesić, and I. Petrović, “Von Mises Mixture PHD Filter,” _IEEE Signal Processing Letters_ , vol. 22, no. 12, pp. 2229–2233, Dec. 2015.
  * [163] Y. Liu, W. Wang, J. Chambers, V. Kilic, and A. Hilton, “Particle Flow SMC-PHD Filter for Audio-Visual Multi-speaker Tracking,” in _Latent Variable Analysis and Signal Separation_ , P. Tichavský, M. Babaie-Zadeh, O. J. Michel, and N. Thirion-Moreau, Eds. Cham: Springer International Publishing, 2017, pp. 344–353.
  * [164] S. Blackman and R. Popoli, _Design and Analysis of Modern Tracking Systems_. Norwood (Massachusetts), USA: Artech House, 1999.
  * [165] H. W. Kuhn, “The Hungarian Method for the Assignment Problem,” _Naval Research Logistics Quarterly_ , vol. 2, pp. 83–97, Mar. 1955.
  * [166] R. L. Rothrock and O. E. Drummond, “Performance Metrics for Multiple-Sensor Multiple-Target Tracking,” _Proc. SPIE_ , vol. 4048, Jul. 2000.
  * [167] A. A. Gorji, R. Tharmarasa, and T. Kirubarajan, “Performance Measures for Multiple Target Tracking Problems,” in _Proc. of Intl. Conf. on Information Fusion_ , Chicago (Illinois), USA, Jul. 2011.
  * [168] D. Schuhmacher, B.-T. Vo, and B.-N. Vo, “A Consistent Metric for Performance Evaluation of Multi-Object Filters,” _IEEE Trans. on Signal Processing_, vol. 56, no. 8, pp. 3447–3457, Aug. 2008.
  * [169] A. S. Rahmathullah, A. F. García-Fernández, and L. Svensson, “Generalized Optimal Sub-Pattern Assignment Metric,” in _International Conference on Information Fusion (Fusion)_ , Jul. 2017, pp. 1–8.
  * [170] B. Ristic, B.-N. Vo, and B.-T. Vo, “A Metric for Performance Evaluation of Multi-Target Tracking Algorithms,” _IEEE Trans. on Signal Processing_, vol. 59, no. 7, pp. 3452–3457, Jul. 2011.
  * [171] T. Gerkmann and R. C. Hendriks, “Unbiased MMSE-Based Noise Power Estimation With Low Complexity and Low Tracking Delay,” _IEEE/ACM Trans. on Audio, Speech, and Language Processing_ , vol. 20, no. 4, pp. 1383–1393, May 2012\.

|  Christine Evers (M’14-SM’16) is a lecturer in Computer Science at the
University of Southampton. She was the recipient of an EPSRC Fellowship,
hosted at Imperial College London between 2017-2019. She worked as a research
associate at Imperial College London between 2014-2017; as a senior systems
engineer at Selex Electronic Systems between 2010-2014; and as a research
fellow at the University of Edinburgh between 2009-2010. She received her PhD
from the University of Edinburgh in 2010; her MSc degree in Signal Processing
and Communications from the University of Edinburgh in 2006; and her BSc
degree in Electrical Engineering and Computer Science from Jacobs University,
Germany, in 2005. Her research focuses on Bayesian learning for machine
listening. She is currently member of the IEEE SPS Technical Committee on
Audio and Acoustic Signal Processing and serves as an associate editor of the
EURASIP Journal on Audio, Speech, and Music Processing.
---|---
|  Heinrich W. Löllmann (M’08-SM’18) received the Dipl.-Ing. (univ) degree in
electrical engineering in 2001 and the Dr.-Ing. degree in 2011 at RWTH Aachen
University. He worked as scientific co-worker at the Institute of
Communication Systems and Data Processing of RWTH Aachen University from 2001
to 2012. Since 2012, he is a senior researcher at the Chair of Multimedia
Communications and Signal Processing at the Friedrich-Alexander University
Erlangen-Nürnberg (FAU). Heinrich Löllmann has authored one book chapter and
more than 40 refereed papers. His research focuses on speech and audio signal
processing, including filter-bank design, speech dereverberation and noise
reduction, estimation of room acoustical parameters, and algorithms for robot
audition. He is currently member of the IEEE SPS Technical Committee on Audio
and Acoustic Signal Processing.
---|---
|  Heinrich Mellmann studied Computer Science and Mathematics at the Humboldt-
Universität zu Berlin (HUB), and received his degrees Dipl.-Inf. in 2011 and
Dipl.-Math. in 2017 respectively. He is currently a research assistant at the
chair of Adaptive Systems, HUB, pursuing a PhD in the area of cognitive
robotics. Over the past years, he was involved in a number of research
projects in cognitive robotics. Heinrich Mellmann has been actively involved
in the RoboCup initiative since 2005, and has been leading the RoboCup team
‘Berlin United’ since 2009. His research interests within cognitive robotics
focus on spatial perception and decision making in autonomous humanoid robots.
---|---
|  Alexander Schmidt (M’19) studied Electrical Engineering at Friedrich-
Alexander-University Erlangen-Nürnberg (FAU), Germany, and received his MSc
degree in 2015. He is currently a research assistant at the Chair of
Multimedia Communications and Signal Processing, FAU, pursuing a PhD in the
area of multichannel signal enhancement for robot audition. He was with the EU
FP7 Project ‘Embodied Audition for Robots (EARS)’. His special interests lie
in the area of (sparse) dictionary learning for signal representation combined
with physical-mechanical models.
---|---
|  Hendrik Barfuss received a Master degree from Friedrich-Alexander-
University Erlangen-Nürnberg (FAU), Germany, in 2013. He is currently pursuing
a Dr.-Ing. degree in Electrical, Electronic and Communication Engineering at
the Chair of Multimedia Communications and Signal Processing, FAU. His
research interests are in microphone array signal processing, spatial
filtering, speech enhancement, and localization.
---|---
|  Patrick A. Naylor (M’89, SM’07, F’20) is Professor of Speech and Acoustic
Signal Processing at Imperial College London. He received the BEng degree in
Electronic and Electrical Engineering from the University of Sheffield, UK,
and the PhD degree from Imperial College London, UK. His research interests
are in speech, audio and acoustic signal processing. His current research
addresses microphone array signal processing, speaker diarization, and
multichannel speech enhancement. He has also worked on speech dereverberation
including blind multichannel system identification and equalization, acoustic
echo control, non-intrusive speech quality estimation, and speech production
modelling with a focus on the analysis of the voice source signal. In addition
to his academic research, he enjoys several collaborative links with industry.
He is currently a member of the Board of Governors of the IEEE Signal
Processing Society and President of the European Association for Signal
Processing (EURASIP). He was formerly Chair of the IEEE Signal Processing
Society Technical Committee on Audio and Acoustic Signal Processing. He has
served as an associate editor of IEEE Signal Processing Letters and is
currently a Senior Area Editor of IEEE Transactions on Audio Speech and
Language Processing.
---|---
|  Walter Kellermann is a professor for communications at the Friedrich-
Alexander-University University Erlangen-Nürnberg (FAU), Germany, since 1999.
He received the Dipl.-Ing. (univ.) degree in Electrical Engineering from FAU,
in 1983, and the Dr.-Ing. degree from the Technical University Darmstadt,
Germany, in 1988. From 1989 to 1990, he was a postdoctoral Member of Technical
Staff at AT&T Bell Laboratories, Murray Hill, NJ. In 1990, he joined Philips
Kommunikations Industrie, Nuremberg, Germany, to work on hands-free
communication in cars. From 1993 to 1999, he was a Professor at the
Fachhochschule Regensburg, where he also became Director of the Institute of
Applied Research in 1997. In 1999, he cofounded DSP Solutions, a consulting
firm in digital signal processing, and he joined FAU as a Professor and Head
of the Audio Research Laboratory. He authored or coauthored 21 book chapters,
300+ refereed papers in journals and conference proceedings, as well as 70+
patents, and is a co-recipient of ten best paper awards. His current research
interests include speech signal processing, array signal processing, adaptive
and learning algorithms, and its applications to acoustic human–machine
interfaces. Dr. Kellermann served as an Associate Editor and Guest Editor to
various journals, including the IEEE Transactions on Speech and Audio
Processing from 2000 to 2004, the IEEE Signal Processing Magazine in 2015, and
presently serves as Associate Editor to the EURASIP Journal on Applied Signal
Processing. He was the General Chair of seven mostly IEEE-sponsored workshops
and conferences. He served as a Distinguished Lecturer of the IEEE Signal
Processing Society (SPS) from 2007 to 2008. He was the Chair of the IEEE SPS
Technical Committee for Audio and Acoustic Signal Processing from 2008 to
2010, a Member of the IEEE James L. Flanagan Award Committee from 2011 to
2014, a Member of the SPS Board of Governors (2013-2015), Vice President
Technical Directions of the IEEE Signal Processing Society (2016-2018) and is
currently a member of the SPS Nominations Appointments Committee (2019-2020).
He was awarded the Julius von Haast Fellowship by the Royal Society of New
Zealand in 2012 and the Group Technical Achievement Award of the European
Association for Signal Processing (EURASIP) in 2015. In 2016, he was a
Visiting Fellow at Australian National University, Canberra, Australia. He is
an IEEE Fellow.
---|---
  *[BSS]: Blind Source Separation
  *[SLAM]: Simultaneous Localization and Mapping
  *[LOCATA]: LOCalization And TrAcking
  *[DoA]: Direction-of-Arrival
  *[DICIT]: Distant talking Interfaces for Control of Interactive TV
  *[CSTR]: Centre for Speech Technology Research
  *[VCTK]: Voice Cloning ToolKit
  *[IR]: Infra-Red
  *[VAP]: Voice-Active Period
  *[MUSIC]: MUltiple SIgnal Classification
  *[TDE]: Time Delay Estimation
  *[TDoA]: Time Delay of Arrival
  *[DTFT]: Discrete-Time Fourier Transform
  *[GCC]: Generalized Cross-Correlation
  *[PHAT]: PHAse Transform
  *[ToA]: Time-of-Arrival
  *[AIR]: Acoustic Impulse Response
  *[EVD]: EigenValue Decomposition
  *[LMS]: Least-Mean-Square
  *[RTF]: Relative Transfer Function
  *[RTFs]: Relative Transfer Function
  *[LS]: Least Squares
  *[EM]: Expectation-Maximization
  *[TDoAs]: Time Delay of Arrival
  *[ITDs]: Interaural Time Difference
  *[ILDs]: Interaural Level Difference
  *[IPDs]: Interaural Phase Difference
  *[SRP]: Steered Response Power
  *[MVDR]: Minimum Variance Distortionless Response
  *[ESPRIT]: Estimation of Signal Parameters via Rotational Invariance Techniques
  *[SVD]: Singular Value Decomposition
  *[PHD]: Probability Hypothesis Density
  *[STFT]: Short-Time Fourier Transform
  *[]:
  *[GEVD]: Generalized   ()
  *[ILD]: Interaural Level Difference
  *[IPD]: Interaural Phase Difference
  *[DNN]: Deep Neural Network
  *[pdf]: Probability Density Function
  *[CDR]: Coherent-to-Diffuse Ratio
  *[PDAF]: Probabilistic Data Association Filter
  *[JPDA]: Joint Probabilistic Data Association
  *[RFS]: Random Finite Set
  *[RFSs]: Random Finite Set
  *[ID]: Identity
  *[VAPs]: Voice-Active Period
  *[FAR]: False Alarm Rate
  *[VAD]: Voice Activity Detector
  *[TL]: Track Latency
  *[TFR]: Track Fragmentation Rate
  *[OSPA]: Optimal SubPattern Assignment
  *[IDs]: Identity
  *[DRR]: Direct-to-Reverberant Ratio
  *[CPSD]: Cross-Power Spectral Density
  *[FARs]: False Alarm Rate
  *[PSD]: Power Spectral Density
  *[DPD]: Direct Path Dominance
  *[TFRs]: Track Fragmentation Rate