Fibroblast mediated dynamics in diffusively uncoupled myocytes - a simulation
study using 2-cell motifs.

S. Sridhar 1,* and Richard.H.Clayton 1

1 Department of Computer Science, University of Sheffield, Sheffield, United
Kingdom

*<EMAIL_ADDRESS>
## Abstract

In healthy hearts myocytes are typically coupled to nearest neighbours through
gap junctions. Under pathological conditions such as fibrosis, or in scar
tissue, or across ablation lines myocytes can uncouple from their neighbours.
Electrical conduction may still occur via fibroblasts that not only couple
proximal myocytes but can also couple otherwise unconnected regions. We
hypothesise that such coupling can alter conduction between myocytes via
introduction of delays or by initiation of premature stimuli that can
potentially result in reentry or conduction blocks. To test this hypothesis we
have developed several $2$-cell motifs and investigated the effect of
fibroblast mediated electrical coupling between uncoupled myocytes. We have
identified various regimes of myocyte behaviour that depend on the strength of
gap-junctional conductance, connection topology, and parameters of the myocyte
and fibroblast models. These motifs are useful in developing a mechanistic
understanding of long-distance coupling on myocyte dynamics and enable the
characterisation of interaction between different features such as myocyte and
fibroblast properties, coupling strengths and pacing period. They are
computationally inexpensive and allow for incorporation of spatial effects
such as conduction velocity. They provide a framework for constructing scar
tissue boundaries and enable linking of cellular level interactions with scar
induced arrhythmia.

## Introduction

The heart is an electro-mechanical pump, whose coordinated mechanical
contraction is enabled by the propagation of synchronized waves of electrical
excitation. In the mammalian heart, myocytes and fibroblasts constitute two of
the most important type of cells, with the larger myocytes being responsible
for cardiac electrical activity and the smaller but more numerous fibroblasts
maintaining the electro-mechanical integrity of the heart [1]. Fibroblasts
play a crucial role in the repair of heart muscles especially post injury or
disease[1, 2]. In aged or diseased hearts the number of fibroblasts may
increase substantially (up to $40\%$[3]) resulting in increased collagen
deposition causing fibrosis.

Furthermore fibroblasts themselves differentiate into much larger
myofibroblasts in injured and diseased hearts. In this paper we have used the
terms fibroblasts and myofibroblasts interchangeably even though we are mostly
referring to myofibroblasts. For a long time fibroblasts were not believed to
influence the electrical conduction in heart muscles and even today the exact
nature of myocyte-fibroblast coupling in vivo is debated [2, 4]. However
several in vitro experiments have reported the existence of gap-junctional
coupling between myocytes and fibroblasts under both physiological and
pathological conditions [5, 6, 7, 8]. Experiments have shown that coupling
between myocytes and fibroblasts can significantly alter the conduction
properties of the tissue [5, 9], excitability of myocytes [10] and their
resting membrane potential [6]. At the level of the tissue and organ, fibrosis
is known to substantially affect wave propagation in the heart and it is well
understood that fibrosis can create a substrate for cardiac arrhythmia [11,
12, 13, 14, 15, 16].

Several in silico studies have investigated how fibroblasts alter the
electrical activity of individual myocytes and tissue under both normal and
pathological conditions [17, 18, 19, 20, 21, 22, 23]. These studies have
usually simulated fibroblasts as being either attached to a myocyte or
inserted in a tissue of coupled myocytes thereby coupling nearest neighbours.
However heterocellular cell culture experiments have shown that fibroblasts
can enable conduction up to $300\mu m$ [5]. In vivo fibroblasts have been
observed to form large sheet-like extensions having additional folds and
elongated cytoplasmic processes [2]. In the sinoatrial node it has been
observed that an individual fibroblast could form membrane juxtaposition with
a nearby myocytes covering up to $720$ $\mu m^{2}$ [2, 24]. Fibroblasts that
have such long extensions can potentially couple with multiple myocytes that
are spatially distant.

Such long range interactions between distant myocytes mediated via fibroblasts
have the potential to modify tissue electrophysiology and dynamics via
conduction delays and by exciting resting and partially recovered regions. The
possibility of such long range coupling increases in diseased or injured
hearts where there are a larger number of myofibroblasts. Ablation lines are
another scenario where the fibroblasts that aid the repair of the ablation
scars could produce conduction pathways between regions that have been
electrically isolated [25]. Another possible scenario occurs when islands of
myocytes are trapped in a scar within a sea of fibroblasts resulting in
conduction pathways that connect distant uncoupled regions [26, 27].

Most computational studies have approached the problem of fibroblast induced
dynamics in terms of the effect of the fibroblast distribution, their density
[22, 23, 28] and texture [29, 30]. On the other hand in this paper we describe
a bottom-up approach by considering several simple motifs consisting of $1$ or
$2$ fibroblast units coupled to a pair of mutually uncoupled myocytes.The idea
of motifs has previously been used successfully to understand networks in
complex biological domains including gene transcription, biochemical systems,
neuronal networks [31, 32, 33]. In this paper we have extended the idea of
motif to model complex interactions between myocytes and fibroblasts at the
scar boundaries.

We hypothesise that fibroblast mediated coupling that connects diffusively
uncoupled myocytes can potentially alter conduction between myocytes via the
introduction of delays or by the initiation of premature stimuli and can
potentially result in reentry or conduction blocks in tissue. We have used the
different motifs to investigate the effect of different configurations of
fibroblast mediated electrical coupling between mutually uncoupled myocytes.
We have also determined the effect of the different motif features such as
connection topology, coupling strength, pacing period and cell parameters on
the dynamics of the individual myocytes.

The $2$-cell motifs allow for easy comparison of the influence of the
intrinsic myocyte dynamics (characterized by its restitution) with the local
coupling topology (characterized by the individual motifs) and pacing period.
The advantage of using motifs is that they allow for a computationally
inexpensive approach to study the effect of individual variation of the many
features of the myocyte-fibroblast coupled system. By incorporating delay in
the stimulation times of individual myocytes of the motif, we study the
combined effect of conduction velocity in the tissue and the conduction delay
arising from the coupling of disconnected myocytes via fibroblasts.
Furthermore we have characterized the regimes of myocyte dynamics that can
arise from this kind of fibroblast mediated long-range conduction.

## Methods

The electrical activity of myocytes was described using the $TNNP-TP06$ model
of human ventricular cells [34, 35], while the electrophysiological properties
of the fibroblasts were described using the $MacCannell$ “active” fibroblast
model [36]. The time variation of the transmembrane voltage $V$ for myocytes
coupled to $n$ fibroblasts was given as,

$C_{m}\times\frac{dV}{dt}=-I_{ion}+\sum_{k=1}^{n}G_{gap}(V-V_{f})$ (1)

Here $I_{ion}$ describes the total of all ionic currents:

$I_{ion}=I_{Na}+I_{to}+I_{K1}+I_{Kr}+I_{Ks}+I_{CaL}+I_{NaCa}+I_{NaK}+I_{pCa}+I_{pK}+I_{bCa}+I_{bNa}$
(2)

where $I_{Na}$ is the sodium current, $I_{to}$ is the transient outward
current, $I_{K1}$, $I_{Kr}$ and $I_{Ks}$ are the inward rectifier, delayed
rectifier and slow delayed rectifier potassium currents, $I_{CaL}$ is the
L-type $Ca^{2+}$ current, $I_{NaK}$ is the $Na^{+}/K^{+}$ pump current,
$I_{NaCa}$ is the $Na^{+}/Ca^{2+}$ exchanger current, $I_{pCa}$ and $I_{pK}$
plateau calcium and potassium currents and $I_{bCa}$ and $I_{bNa}$ are the
background $Na^{+}$ and $Ca^{2+}$ currents. $V_{f}$ is the fibroblast
transmembrane potential while $G_{gap}$ is the strength of the gap junctional
coupling between myocyte and fibroblast.

The $MacCannell$ fibroblast model equations are used to describe the time
evolution of the fibroblast membrane potential $V_{f}$ (similar to Equation. 1
and Equation. 2) with the ionic currents comprised of inward rectifying
potassium current $I_{fK1}$, the time- and voltage -dependent potassium
currents $I_{fKv}$, $I_{fNaK}$ a sodium-potassium pump current and a
background sodium current $I_{bNa^{+}}$. $C_{m}$ and $C_{f}$ are the cell
capacitance per unit surface area of myocyte and fibroblast set to $150$pF and
$50$pF (corresponding to the larger myofibroblast [37]) respectively.

For the myocytes, we used two parameter sets corresponding to Shallow and
Steep restitution slopes (see Table $2$, slope $=0.7$ for Shallow and slope
$=1.8$ for Steep in ten Tusscher et al [35]). The uncoupled fibroblast resting
membrane potential $V_{FR}$ were set to either $-24.5$ mV or $-49.0$ mV [12].
Most of the results described in the paper were obtained with $V_{FR}$ set to
$-24.5$ mV. A subset of results with $V_{FR}=-49.0$ mV is described in the
supplementary material. The different resting membrane potentials were
obtained by shifting the gating variable voltage dependence of the time
dependent potassium current [19].

The action potential of the uncoupled myocyte stimulated at a period of $600$
ms for both Shallow and Steep parameters is shown in Fig. 1(a) while Fig.
1(b,d) show the effect of weak fibroblast coupling on the myocyte and
fibroblast trans-membrane potentials respectively. In Fig. 1(c) we have
plotted the restitution curve for the Shallow and Steep uncoupled myocytes.

Figure 1: Time series of transmembrane potentials for myocytes and
fibroblasts. (a) Transmembrane voltage for uncoupled myocyte (a) and the
corresponding S1S2 restitution curves generated at $T=600$ ms pacing for both
Shallow (broken) and Steep (solid) parameters (c). The time series of the
transmembrane potential of myocyte (b) and fibroblast (resting potential
$V_{FR}=-49$ mV) (d) for the case of weak coupling ($G_{gap}=0.5$ nS) for both
Shallow (broken) and Steep (solid) parameters.

In our study we considered two identical myocytes coupled either to $1$ or $2$
fibroblast units via different connection motifs (Fig. 2). Each fibroblast
unit represents a fixed number of fibroblasts (see below) connected in
parallel and coupled to a neighbouring myocyte (similar to Figure.$1$ in [36]
and Figure.$2$ in [38]). The different topology considered are (i) $Motif-1$
where one fibroblast unit is coupled to two myocytes, (ii) $Motif-2$ where
there is an asymmetry in the number of passive cells connected to each myocyte
and (iii) $Motif-3$ where two fibroblast units are coupled in a symmetric
fashion with the two myocytes. In the following text we have referred to each
fibroblast unit as a single fibroblast (viz., Fib $1$ and Fib $2$ in Fig. 2).
The myocyte-fibroblast interaction is modelled as occurring via a gap-
junction. We considered two kinds of gap-junction conductance viz., $G_{Loc}$
and $G_{Long}$ corresponding to electrical conduction between fibroblast and
the proximal and distal myocytes respectively. Note that for the results
reported here the coupling strength of a fibroblast to a proximal myocyte
$G_{Loc}$, if a coupling exists between them, is the same for both fibroblast
units. Similarly if coupling exists the coupling strength of a fibroblast to a
distal myocyte $G_{Long}$ is the same for both myocyte fibroblast units. For
the gap junctional coupling we have chosen values falling in the range ($0-4$
nS) that is considered to be representative of the effect of fibroblasts in
cell-cultures [19].

We captured the effect of conduction velocity in tissue by the introduction of
a delay ($\tau_{D}$) in the stimulation of the distal myocyte (i.e., the
myocyte coupled to a fibroblast with a coupling strength of $G_{Long}$). While
propagation delays of $11-68$ ms have been observed in heterocellular culture
[5], in our study we performed simulations with two representative time delays
of $\tau_{D}=10$ ms and $\tau_{D}=25$ ms.

Figure 2: Connection topology. The different motifs constructed with two
myocytes coupled gap-junctionally to two fibroblast units, viz., $Motif-1$
(a), $Motif-2$ (b) and $Motif-3$ (c).

The resulting equation for the coupled system $Motif-3$ (Fig. 2 (c)) is:

$\frac{dV_{i}}{dt}=-I_{ion_{i}}+\frac{N_{f}}{C_{m}}\times(G_{Loc_{i}}\times(V_{f_{i}}-V_{i})+G_{Long_{i}}\times(V_{f_{j}}-V_{i}))$
(3)
$\frac{dV_{f_{i}}}{dt}=-I_{fibion_{i}}+\frac{1}{C_{f}}\times(G_{Loc_{i}}\times(V_{i}-V_{f_{i}})+G_{Long_{j}}\times(V_{j}-V_{f_{i}}))$
(4)

where $i={1,2}$, $j={1,2}$ and $i\neq j$. Here $V_{i}$ and $V_{f_{i}}$
represent the transmembrane voltage of the $i^{th}$ myocyte and $i^{th}$
fibroblast respectively, while $N_{f}$ represents the fibroblast density
(number of fibroblasts in a fibroblast unit). The equations for $Motif-1$ are
obtained by setting $G_{Loc1}$ and $G_{Long2}$ to zero. Similarly $Motif-2$
equations are obtained by setting $G_{Long2}=0$ nS.

The coupled ordinary differential equations were solved using an explicit
adaptive time step method with a maximum time-step of $0.001$ ms [39]. In
order to obtain an action potential, either or both the myocytes were
individually stimulated every $T$ ms by applying an external stimulus of
strength $-52~{}\mu A/mm^{2}$ for a duration of $1$ ms. While many of the
results described in the paper were obtained for $T=600$ ms, the effect of
pacing period was also investigated by setting $T=500,400$ and $300$ ms
respectively. The fibroblast density $N_{f}$ was set to $4$. Each simulation
was performed by stimulating the myocytes $20$ times, the first $10$ action
potentials generated were ignored and the mean action potential duration and
fraction of action potentials in the Pacing and Response cells were calculated
over the last $10$ action potentials alone. The initial conditions of the
variables describing both myocytes and fibroblasts were set to their uncoupled
resting membrane values.

## Results

### Effect of fibroblast coupling on APD

Figure 3: Effect of coupling on APD. Two-parameter plots for all the $3$
motifs describing the effect of the coupling strengths $G_{Loc},G_{Long}$ on
the APD of both $myocyte-1$ (a, c, e) and $myocyte-2$ (b,d,f) with Steep
parameters.

In order to characterize the effect of connection topology and strength of
myocyte-fibroblast coupling, we plotted $2$-parameter portraits for the mean
APD of the myocytes. The two parameters $G_{Loc}$ and $G_{Long}$ were varied
over the range $0-4$ nS in steps of $0.5$ nS. Figure. 3 shows the effect of
strength of coupling on mean APD for both $myocyte-1$ (Fig. 3 (a,c,e) and
$myocyte-2$ (Fig. 3 (b,d,f)) for the different motifs with Steep restitution
parameters. (See Supplementary Fig. 10 for the equivalent figure for Shallow
parameters).

While the general effect of fibroblast coupling was to reduce the myocyte
$APD$, there was significant variation across motifs. $Motif-1$ was the
simplest motif with one fibroblast coupled to two myocytes, and showed a
reduction in APD of nearly $70$ ms for the case of the strongest coupling
($G_{Loc}=G_{Long}=4.0$ nS) as compared to the case of no coupling (Fig. 3
(a,b)). However the largest reduction in APD in $myocyte-1$ occurred for the
case of ($G_{Loc},G_{Long}$) $=(0.0,4.0)$ nS (For $myocyte-2$,
($G_{Loc},G_{Long}$) $=(4.0,0.0)$ nS). The reduction in APD for Shallow
parameters was smaller than that of the Steep parameters, with a maximum
reduction of APD compared to that of no coupling being around $44$ ms (see
Supplementary Fig. 10). In the case of Shallow parameters the maximum
reduction occurred for the case of ($G_{Loc}=G_{Long}=4.0$ nS).

We next considered the effect of APD for $Motif-2$ which has a structural
asymmetry with each myocyte exposed to a different number of fibroblasts.
Figure. 3(c,d) shows the effect of $Motif-2$ coupling on APD for both
$myocyte-1$ and $myocyte-2$ for the Steep parameters. Due to the asymmetry in
their coupling the magnitude of APD was not the same for both myocytes for a
given ($G_{Loc},G_{Long}$) pair. However APD of both myocytes became equal as
$G_{Long}$ decreased. The reduction of APD in $myocyte-1$ was greater than in
$myocyte-2$. This motif too showed a greater reduction in APD for the Steep
parameters as compared to the Shallow parameters (Supplementary Fig. 10(c,d)).
$Motif-3$ is symmetric in terms of coupling and the reduction in APD for a
given coupling strength was the same for both myocytes, i.e. the difference in
APD between the myocytes was zero. However the decrease in APD with coupling
strength was larger for the Steep parameter set (Fig. 3 (e,f)) compared to the
Shallow parameter case(Supplementary Fig. 10 (e,f)). Comparing across the
different motifs we find that the reduction in APD with coupling increased
with complexity of the topology, with the largest reduction in APD occurring
for the case of $Motif-3$.

For a representative set of coupling strengths we have compared the action
potential profiles of both $myocyte-1$ and $myocyte-2$ in $Motif-2$ for both
Shallow and Steep parameters in Fig. 4. For each of the motifs, there was
significant change in features of the myocyte action potential shape including
peak, dome and duration of recovery depending on the coupling strengths
$(G_{Loc},G_{Long})$. The corresponding action potential profiles for
$Motif-1$ and $Motif-3$ are plotted in the Supplementary section (Fig. 11 and
Fig. 12).

Figure 4: Action potential time series for coupled myocytes. Action potential
profiles for both myocytes in $Motif-2$ while pacing the cells at $T=600$ ms
at $G_{Loc}=G_{Long}=0.5$ nS (a), $2.0$ nS (b) and $4.0$ nS (c) respectively.
The solid and broken blue lines correspond to $myocyte-1$ and $myocyte-2$ for
the Shallow parameter set while solid and broken red traces correspond to the
same for the Steep parameter set.

### Effect of delay $\tau_{D}$ in stimulation

Figure 5: Effect of time delay in the stimulation of $myocyte-2$. The
difference in APD ($\Delta APD$) between the two myocytes for $Motif-1$ (a, b,
c), $Motif-2$ (d,e,f) and $Motif-3$ (g, h, i) for Steep parameters. The values
of the delay are $\tau_{D}=0$ ms (a,d,g), $=10$ ms (b, e, h) and $=25$ ms (c,
f, i).

In the result described above, both myocytes were stimulated simultaneously.
To capture the effect of conduction velocity across tissue we stimulated
$myocyte-2$ with a delay after the stimulation of $myocyte-1$. Supplementary
Fig. 14 describes the time series of myocyte membrane voltage with $myocyte-2$
stimulated with a time delay $\tau_{D}=25$ ms. Figure. 5 shows the effect of
stimulating $myocyte-2$ with a delay of $0$ ms, $10$ ms and $25$ ms after
stimulating $myocyte-1$ on $\Delta APD$, (defined as the difference in APD
between $myocyte-1$ and $myocyte-2$). With $Motif-1$, for all values of delay
considered here, we observed that for both Steep and Shallow parameter sets
(Supplementary figure Fig. 13(a)) $\Delta APD$ was positive (negative) for
$G_{Long}=0$ nS ($G_{Loc}=0$ nS). However for the Steep parameters the
magnitude of $\Delta APD$ became more negative or less positive (less negative
or more positive) with increase of $G_{Long}$ ($G_{Loc}$). On the other hand
for the Shallow parameter set the change in the magnitude of $\Delta APD$ was
not monotonic for all values of delay.

We next considered the effect of delay for the case of $Motif-2$. In the case
of zero delay in stimulation for Shallow parameters $myocyte-2$ showed a
greater APD than $myocyte-1$ only for $G_{Loc}\leq 1nS$. For stronger local
coupling ($G_{Loc}>1$ nS), $\Delta APD$ was positive indicating that the
reduction in APD of $myocyte-2$ due to coupling was greater than the reduction
in that of $myocyte-1$ (Supplementary figure Fig. 13 (d)). However for the
case of Steep parameters with no delay, $\Delta APD\geq 0$ only for
$G_{Long}=0$ or when $G_{Loc}$ was large (Fig. 5(d)). For small $G_{Loc}$ the
reduction in APD was greater for $myocyte-1$. Furthermore the magnitude of
$\Delta APD$ was more negative for Steep compared to the Shallow parameter
set. While for stronger coupling, the magnitude of $\Delta APD$ became more
positive for longer delays for both Shallow and Steep parameters, with a
greater increase for Shallow compared to Steep parameters (Supplementary Fig.
13 (e,f)). On the other hand for the Steep case (Fig. 5 (e,f), the reduction
in APD of $myocyte-1$ compared to $myocyte-2$ for small $G_{Loc}$ increased
with the magnitude of delay. However for all values of delay, $\Delta APD<0$
ms occurred only when $G_{Loc}\leq 1.5$ nS.

In $Motif-3$, for a delay of $0$ ms both the myocytes had the same $APD$ as
seen in Fig. 5 (g). For a delay of $10$ ms, both Steep (Fig. 5(h)) and Shallow
parameters (Supplementary figure Fig. 13(h)) show a reduction in APD for
$myocyte-2$ (corresponding to a positive value of $\Delta APD$) for all pairs
of ($G_{Loc},G_{Long}$). However for a delay of $25$ ms alone the Shallow
parameter set showed a reduction in $APD$ in $myocyte-2$ for all conductance
pairs (Supplementary Fig. 13(i)). On the other hand for the Steep parameter
set (Fig. 5(i)) $\Delta APD$ was positive at large ($G_{Loc},G_{Long}$)
values, while a large $G_{Long}$ and small $G_{Loc}$ conductance results in a
greater reduction of APD for $myocyte-1$ than for $myocyte-2$ (as seen by the
negative $\Delta APD$ values).

### Initiating action potentials in a resting cell via fibroblasts

We next considered the case of infinite delay in stimulating one of the cells.
In other words only one myocyte in a motif was stimulated (Pacing Cell) and
the effect of this stimulus on the other myocyte (Response Cell) via the
fibroblast was determined. The goal was to identify regions in the
$2$-parameter conductance space that would result in excitation of the
quiescent Response Cell. In this experiment, the Pacing Cell was stimulated at
a fixed period and the effect on the Response Cell was characterized in terms
of the fraction of stimuli in the Pacing Cell that elicited a response in the
Response Cell.

In Figure. 6 (a-b), the black region $(NR)$ corresponds to conductance value
(small $G_{Loc}$ or $G_{Long}$) that did not produce any excitation in the
Response Cell, even though the Pacing Cell produced an action potential for
every applied stimulus. However for sufficiently strong coupling (larger
values of $G_{Loc}$ and $G_{Long}$) every excitation in the Pacing Cell, was
followed by an excitation in the Response Cell. This parameter region was
marked $1:1$. Between $NR$ and $1:1$ regions, there was an intermediate region
of parameter values where the Response Cell did not respond to every
excitation of Pacing Cell. Rather there was an intermediate response $(IR)$,
with only some of the action potentials in Pacing Cell resulting in an action
potential in the Response Cell.

Figure. 6(c-f) shows the action potentials for the set of parameter points
marked $c,d,e,f$ on Fig. 6(a). The initiation of excitation at a distal (not
directly coupled) myocyte was critically dependent on the ($G_{Loc},G_{Long}$)
conductance values in the motif. For $G_{Long}=0.5$ nS there was no response
(Fig. 6)(c). When the value of $G_{Long}$ was further increased to $1.0$, a
shorter action potential was elicited in Response Cell for every stimulation
of the Pacing Cell (Fig. 6)(d). However further increase of $G_{Long}$ to
$1.5$ nS gave rise to more complex intermittent dynamics in the Response Cell
(Fig. 6)(e). Finally for $G_{Long}\geq 2.5$ nS, the Pacing Cell elicited a
$1:1$ response in the resting cell (Fig. 6)(f).

For the cases where there was a depolarization in the Response Cell, there was
always a time delay with respect to the depolarization of Pacing Cell. However
this delay was a function of the conductance strength (as seen in (Fig.
6)(d,f). The action potential profiles in (Fig. 6)(d,f) suggest a reduction in
APD of the Response Cell compared to the Pacing Cell. Furthermore for large
conductance there was an overlap in the repolarization time-series of both
myocytes suggesting a synchronization of the recovery process in both the
myocytes. Similar $2$-parameter plots for the Pacing and Response cells are
plotted for both $Motif-1$ and $Motif-2$ (Supplementary figures Fig. 15 and
Fig. 16).

Since $Motif-2$ had an asymmetry, we considered two scenarios, viz.,
stimulating either $myocyte-2$ (Supplementary figures Fig. 15(c,d) and Fig. 16
(c,d)) or $myocyte-1$ (Supplementary figures Fig. 15(e-f) and Fig. 16(e-f)).
Stimulating only $myocyte-2$ could initiate action potentials in $myocyte-1$
for a large range of conductance for both Steep and Shallow parameter sets.
While for both Steep and Shallow parameters, small $G_{Loc}$ or $G_{Long}$
values did not produce a response in $myocyte-1$, at larger conductance $1:1$
response was obtained for both the parameter sets. The region of intermediate
response $(IR)$ was larger for Steep compared to Shallow. On the other hand a
very small number of conductance pairs could initiate action potential in
$myocyte-2$ when $myocyte-1$ alone was stimulated for the Shallow parameter
set. For the Steep parameter set, $1:1$ response was not observed for any of
the conductance values.

Figure 6: Dynamical regimes described by two-parameter conductance maps.
Dynamical regimes characterised as the fraction of stimulus that elicit an
action potential in the Response Cell for both Shallow (a) and Steep (b) cell
parameters for $Motif-3$. (c-f) The action potential profiles for the
parameters marked (c,d,e,f) in the case of Shallow parameter set (a) comparing
the different dynamical behaviour in the Pacing and Response cells.

In order to identify the effect of the pacing period on the initiation of
action potentials via fibroblasts in a resting cell, we performed the above
simulations for pacing periods $T=500$ ms, $T=400$ ms and $T=300$ ms. The
parameter space diagrams for $T=500$ ms and $T=400$ ms (not shown) were found
to be qualitatively similar to those obtained for $T=600$ ms, with only the
boundaries between the different regimes and size of regimes varying
marginally depending on the pacing period. On the other hand pacing at $T=300$
ms produced regimes that were spatially more patchy suggesting sensitive
dependence to coupling strength at very rapid pacing. Furthermore at $T=300$
ms pacing $IR$ regimes were also observed in the Pacing Cell implying that at
very rapid pacing $1:1$ response was not always guaranteed especially when
coupled to other cells that can act as a current sink. For $T=300$,
supplementary figures Fig. 17 and Fig. 18 describe the different dynamical
regimes obtained by stimulating either $myocyte-1$ or $myocyte-2$ for the case
of $Motif-2$. Supplementary Fig. 19 and Fig. 20 describe the different regimes
obtained by pacing one cell with $T=300$ ms for $Motif-1$ and $Motif-3$.

Figure 7: Effect of restitution. Fraction of occurrence of each dynamical
regime, viz., $NR$, $IR$ and $1:1$ summed over all pacing cyles and motifs at
both Pacing Cell ($PCell$) and Response Cell ($RCell$) for Shallow (a) and
Steep (b) parameter sets.

The $2$-parameter plots are useful in identifying the
different dynamical regimes and their relation to myocyte-fibroblast coupling
strength. However in order to determine the influence of each of the features
(viz., connection topology, myocyte parameters, pacing period) on the myocyte
dynamics, we determined the fraction of instances for every regime, keeping
one feature fixed at a time. In Fig. 7 we plot the fraction of occurrence of
each regime for Steep and Shallow parameters individually, summing across the
other features for both Pacing Cell ($PCell$) and Response Cell ($RCell$).
Similarly Fig. 8 and Fig. 9 describe the fraction of occurrence of the
different regimes for the different pacing periods and motifs respectively.

Figure 8: Effect of pacing period. Fraction of occurrence of each dynamical
regime, viz., $NR$, $IR$ and $1:1$ summed over all motifs and restitution
types at both Pacing Cell ($PCell$) and Response Cell ($RCell$) for pacing
periods $T=600$ ms, $T=500$ ms, $T=400$ ms and $T=300$ ms.

Based on Fig. 7, Fig. 8 and Fig. 9, we can summarise that (i) Obtaining $1:1$
response was more difficult with Steep parameters compared to Shallow. On the
other hand Steep showed $IR$ for more parameters than Shallow. (ii) With
increase in pacing period, there were fewer instances of $1:1$ and more cases
of $IR$. For $T\leq 400$ ms the Pacing Cell can also elicit $IR$ instead of
$1:1$ for some coupling strengths. (iii) $Motif-3$ with the smallest region of
$NR$ and most instances of $1:1$ allowed for conduction over the largest range
of parameters. On the other hand $Motif-2$ ($T2=0$ in Fig. 2(b)) with the
maximum instances of $NR$ had the smallest fraction of parameters that allowed
any conduction ($IR$ or $1:1$). $Motif-1$ had the least number of $1:1$ among
all motifs but had a larger fraction of cases describing $IR$ than $Motif-2$.

Figure 9: Effect of connection topology. Fraction of occurrence of each
dynamical regime, viz., $NR$, $IR$ and $1:1$ summed over all restitution types
and pacing cycle lengths at both Pacing Cell ($PCell$) and Response Cell
($RCell$) for $Motif-1$ (a), $Motif-2$ (b,c) and $Motif-3$ (d).

We also repeated our simulations with a more negative fibroblast resting
membrane potential $V_{FR}=-49.0$ mV for $Motif-2$ and $Motif-3$. The results
for a subset of these simulations are described in the Supplementary section.

A more negative fibroblast resting membrane potential resulted in a greater
decrease in APD for both Shallow (Supplementary Fig. 23) and Steep
(Supplementary Fig. 21) restitution parameters compared to a resting membrane
potential of $-24.5$ mV. Supplementary figures Fig. 22 and Fig. 24 describe
the $2$-parameter plots identifying the different dynamical regimes when only
one myocyte was stimulated for Steep and Shallow regimes respectively. While
the different dynamical regimes are qualitatively similar to those obtained in
$V_{FR}=-24.5$mV, it can be observed that due to the more negative fibroblast
resting potential a stronger coupling was required to elicit an action
potential in the Response Cell.

## Discussion

In this paper we have described a simple way to capture complex cellular
dynamics that can occur due to the interaction between $2$ myocytes that are
coupled only via fibroblasts. In the heart the absence of gap-junctional
coupling between myocytes could be due to the presence of ablation lines,
scars or fibrosis resulting in spatial and electrical separation of myocytes.
However under such conditions there is some evidence that conduction can still
occur via fibroblasts across novel pathways coupling otherwise disconnected
myocytes [26, 27, 25, 40]. As we have described in this paper, such
connections can give rise to a range of dynamical behaviour including reduced
APD, synchronization of repolarization in uncoupled myocytes, initiation of
action potentials in resting cells and conduction delays. At the level of
cardiac tissue such non-local coupling might possibly impact wave propagation
and lead to reentry or conduction blocks.

We have considered $3$ different topological arrangements of $2$ myocytes
coupled to either $1$ or $2$ fibroblast units and investigated their effect on
myocyte dynamics. We have chosen myocyte parameters corresponding to Shallow
and Steep restitution slopes [35] and for the fibroblast we have used the
MacCannell “active” fibroblast model [36] with modifications made to obtain
different resting membrane potentials [19]. Since the setup we have considered
here is more likely in a heart tissue undergoing repair from either surgery,
injury or disease we have modified the MacCannell model to simulate
myofibroblasts by setting $C_{f}=50$ pF.

Irrespective of the connection topology, the primary effect of the coupling
was a decrease in the APD of the myocyte with increase in coupling strength
(Fig. 3) with the myofibroblast acting as a leaky capacitor [12]. However the
magnitude of the change in APD depended on the type of motifs and the myocyte
parameters (i.e., whether the restitution is Steep or Shallow). The decrease
in APD is more pronounced for the case with Steep parameters (Fig. 3) compared
to Shallow parameters (Supplementary Fig. 10) across all the motifs. Also the
magnitude of the change in APD increased with the complexity of the
connections with $Motif-3$ (the most complex and tightly coupled of the motifs
considered) showing the greatest decrease in APD at maximum coupling.
Furthermore we observed that the decrease in APD also increases for a more
negative resting membrane potential of the fibroblast (Supplementary figures
Fig. 23 and Fig. 21), an observation consistent with earlier studies [12].

Since the two myocytes considered in the motifs are uncoupled and spatially
separated, there is a time difference in the excitation of the two cells. This
could be due to the difference in times that the cells are excited in the
myocardium (i.e., the effect of conduction velocity across the heart tissue)
and/or due to the effect of conduction delay in the propagation of current via
fibroblasts. We first studied the case where the predominant delay is the time
difference ($\tau_{D}$) between the stimulation of both myocytes (Fig. 5 and
Supplementary Fig. 13). The delay in the stimulation of $myocyte-2$ breaks the
symmetry of coupling and results in differential myocyte APD even for
$Motif-1$ and $Motif-3$ (Supplementary Fig. 14) for the action potential
profiles for cells stimulated at different times). We observed that the sign
and magnitude of the difference in APD of the two myocytes ($\Delta APD$) is
sensitive to one or more of the features varied here viz., connection
topology, coupling strengths, delay ($\tau_{D}$) and the myocyte parameters
(Shallow or Steep). For conditions of finite stimulus delay, $\tau_{D}$ acts
as a proxy to the effect of conduction velocity across the myocardium.

Another type of delay that can occur in these systems is the time taken by the
excitation to propagate from one myocyte to another purely via a fibroblast.
In order to investigate the effect of this conduction we considered the case
of $\tau_{D}\rightarrow\infty$. In practical terms this means that one of the
cells in the motif is never stimulated externally for the duration of the
simulation, while the other myocyte is stimulated periodically as before.
Figure. 6 (a-b) shows the $2$-parameter phase space characterizing the
different dynamical regimes as a function of the gap-junctional coupling
conductance. Figure. 6 (c-f) highlights the sensitive dependence of the
myocyte dynamics on the coupling strength, with small changes in coupling
parameters resulting in widely different myocyte responses. The idea of
infinite delay in stimulation is especially useful to illustrate the dynamics
of initiation of action potential in regions that are isolated from its
neighbours except for conduction via fibroblasts.

The intermittent sequence of action potentials observed in the $IR$ regime in
the Response Cell (and in Pacing Cell at very rapid pacing) suggests a
plausible dynamical mechanism that can result in conduction block and
initiation of reentry in tissue. Fig. 7, Fig. 8 and Fig. 9 summarize the role
of the individual features in determining the dynamical regimes. While steep
restitution and rapid pacing favour complex or irregular dynamics in both the
Response and Pacing Cells, the motif structure and the location of the
stimulated myocyte (with respect to fibroblasts connected to it) also play a
critical role in determining the dynamical regime. In particular $Motif-1$
(where one fibroblast couples to two uncoupled myocytes) with its very low
number of $1:1$ response is a very plausible connection topology in tissue
that can give rise to conduction blocks and reentry.

While heterocellular coupling between myofibroblasts and myocytes have been
reported to initiate ectopic activity in vitro [7], to the best of our
knowledge this is the first in silico study to systematically investigate the
different features that can potentially influence the myocyte dynamics when
coupled purely via fibroblasts. We hypothesised that non-local coupling could
potentially result in the initiation of action potentials in the quiescent
myocyte and identified the parameters that described the various dynamical
regimes possible for different connection topology.

Motifs are a simple prototype to investigate the effect of long-distance
connections between uncoupled myocytes connected only via fibroblasts. Many
studies have looked at the effect of fibroblast distributions in simulated
tissue on wave dynamics [22, 23, 28]. The idea of fibrotic and functional
clusters [41, 42] have been developed based on percolation theory to
investigate the interaction of wave-propagation with fibrosis. Interstitial
fibrosis is associated with non-ischemic cardiomyopathy and has been modelled
as infinitesimal splits in a finite element mesh [43]. Machine learning
algorithms have been employed to understand the effect of local fibrosis
patterns especially at border zones [44]. More recently homogenisation
techniques have been applied to model fibrosis as spatially repeating
structures [45], or by using graph theoretical [46] and volume averaging
approaches [47] to incorporate microscopic structures into a macro-scale
problem.

Our methodology in this study differs significantly from the methods used in
the above papers. We have adopted a bottom-up approach where we have developed
simple structural motifs that can in principle be scaled to build scar
boundary zones. While motifs have been used extensively in other areas [31,
32, 33], to the best of our knowledge this is the first paper that uses the
idea of motifs to describe the various dynamics arising from myocyte-
fibroblast interaction. Our approach is especially suitable to study the
effect of fibroblast connections that form across ablation lines. Fibroblasts
that are involved in the repair of the surgically ablated zones could enable
conduction across the separated regions [25]. The motifs described here
provide a possible structural mechanism to couple disconnected regions in
tissue.

We conclude by stating the limitations of our study and the scope of future
work. While the motifs developed here are a useful prototype to simulate non-
local coupling in cells across ablation line or fibrotic regions, they do not
account for the effect of electrotonic diffusion, which acts to smooth
wavefronts and reduce the effect of local variations in topology. So in as
much as motifs can be a useful tool to quickly explore dynamics while varying
several factors their maximal utility would be realised when they are used to
build scar tissue boundaries and simulate conduction across non-conducting
ablation lines in heart tissue. While we have considered homogeneous myocytes
in our motifs, myocyte properties are heterogeneous in diseased hearts.
Expanding the motif prototype to $2D$ tissue and incorporating heterogeneous
myocytes in the motifs are two areas of future work. In our simulations we
fixed the number of fibroblasts in a unit ($N_{f}$) connected to a myocyte to
be $4$. But the number of fibroblasts coupled to a myocyte is an important
factor that affects the action potential and can be systematically varied to
investigate its effect on wave propagation [48]. Lastly an important
limitation of this paper is that we have only simulated the electrophysiology
of myocytes and not the mechanics of their contraction and relaxation.
Mechanical contractions and the resultant change in tissue geometry have
significant effect on wave-propagation. Another factor that influences wave-
propagation is the mechano-electric feedback. These are aspects that will be
incorporated in future studies.

## References

  * [1] Patrizia Camelliti, Thomas K Borg, and Peter Kohl.  Structural and functional characterisation of cardiac fibroblasts.  Cardiovascular Research, 65(1):40–51, 2005.
  * [2] Peter Kohl and Robert G Gourdie.  Fibroblast–myocyte electrotonic coupling: does it occur in native cardiac tissue?  Journal of Molecular and Cellular Cardiology, 70:37–46, 2014.
  * [3] Tokuhiro Kawara, Richard Derksen, Joris R de Groot, Ruben Coronel, Sara Tasseron, Andre C Linnenbank, Richard NW Hauer, Hans Kirkels, Michiel J Janse, and Jacques MT de Bakker.  Activation delay after premature stimulation in chronically diseased human myocardium relates to the architecture of interstitial fibrosis.  Circulation, 104(25):3069–3075, 2001.
  * [4] Carolina Vasquez, Najate Benamer, and Gregory E Morley.  The cardiac fibroblast: functional and electrophysiological considerations in healthy and diseased hearts.  Journal of Cardiovascular Pharmacology, 57(4):380, 2011.
  * [5] Giedrius Gaudesius, Michele Miragoli, Stuart P Thomas, and Stephan Rohr.  Coupling of cardiac electrical activity over extended distances by fibroblasts of cardiac origin.  Circulation Research, 93(5):421–428, 2003.
  * [6] Michele Miragoli, Giedrius Gaudesius, and Stephan Rohr.  Electrotonic modulation of cardiac impulse conduction by myofibroblasts.  Circulation Research, 98(6):801–810, 2006.
  * [7] Michele Miragoli, Nicolo Salvarani, and Stephan Rohr.  Myofibroblasts induce ectopic activity in cardiac tissue.  Circulation Research, 101(8):755–758, 2007.
  * [8] Lisa Chilton, Wayne R Giles, and Godfrey L Smith.  Evidence of intercellular coupling between co-cultured adult rabbit ventricular myocytes and myofibroblasts.  The Journal of Physiology, 583(1):225–236, 2007.
  * [9] Sharon Zlochiver, Viviana Munoz, Karen L Vikstrom, Steven M Taffet, Omer Berenfeld, and José Jalife.  Electrotonic myofibroblast-to-myocyte coupling increases propensity to reentrant arrhythmias in two-dimensional cardiac monolayers.  Biophysical Journal, 95(9):4469–4480, 2008.
  * [10] E Kizana, SL Ginn, CM Smyth, A Boyd, SP Thomas, DG Allen, DL Ross, and IE Alexander.  Fibroblasts modulate cardiomyocyte excitability: implications for cardiac gene therapy.  Gene Therapy, 13(22):1611–1615, 2006.
  * [11] Kazuhiko Tanaka, Sharon Zlochiver, Karen L Vikstrom, Masatoshi Yamazaki, Javier Moreno, Matthew Klos, Alexey V Zaitsev, Ravi Vaidyanathan, David S Auerbach, Steve Landas, et al.  Spatial distribution of fibrosis governs fibrillation wave dynamics in the posterior left atrium during heart failure.  Circulation Research, 101(8):839–847, 2007.
  * [12] Thao P Nguyen, Yuanfang Xie, Alan Garfinkel, Zhilin Qu, and James N Weiss.  Arrhythmogenic consequences of myofibroblast–myocyte coupling.  Cardiovascular Research, 93(2):242–251, 2012.
  * [13] Norishige Morita, William J Mandel, Yoshinori Kobayashi, and Hrayr S Karagueuzian.  Cardiac fibrosis as a determinant of ventricular tachyarrhythmias.  Journal of Arrhythmia, 30(6):389–394, 2014.
  * [14] Thao P Nguyen, Zhilin Qu, and James N Weiss.  Cardiac fibrosis and arrhythmogenesis: the road to repair is paved with perils.  Journal of Molecular and Cellular Cardiology, 70:83–91, 2014.
  * [15] Gabriel Balaban, Brian P Halliday, Caroline Mendonca Costa, Wenjia Bai, Bradley Porter, Christopher A Rinaldi, Gernot Plank, Daniel Rueckert, Sanjay K Prasad, and Martin J Bishop.  Fibrosis microstructure modulates reentry in non-ischemic dilated cardiomyopathy: insights from imaged guided 2d computational modeling.  Frontiers in Physiology, 9:1832, 2018.
  * [16] Fernando O Campos, John Whitaker, Radhouene Neji, Sébastien Roujol, Mark O’Neill, Gernot Plank, and Martin J Bishop.  Factors promoting conduction slowing as substrates for block and reentry in infarcted hearts.  Biophysical Journal, 117(12):2361–2374, 2019.
  * [17] Vincent Jacquemet.  Pacemaker activity resulting from the coupling with nonexcitable cells.  Physical Review E, 74(1):011908, 2006.
  * [18] Frank B Sachse, Alonso P Moreno, and JA Abildskov.  Electrophysiological modeling of fibroblasts and their interaction with myocytes.  Annals of Biomedical Engineering, 36:41–56, 2008.
  * [19] Vincent Jacquemet and Craig S Henriquez.  Loading effect of fibroblast-myocyte coupling on resting potential, impulse propagation, and repolarization: insights from a microstructure model.  American Journal of Physiology-Heart and Circulatory Physiology, 294(5):H2040–H2052, 2008.
  * [20] Mary M Maleckar, Joseph L Greenstein, Wayne R Giles, and Natalia A Trayanova.  Electrotonic coupling between human atrial myocytes and fibroblasts alters myocyte excitability and repolarization.  Biophysical Journal, 97(8):2179–2190, 2009.
  * [21] Yuanfang Xie, Alan Garfinkel, James N Weiss, and Zhilin Qu.  Cardiac alternans induced by fibroblast-myocyte coupling: mechanistic insights from computational models.  American Journal of Physiology-Heart and Circulatory Physiology, 297(2):H775–H784, 2009.
  * [22] Ivan V Kazbanov, Kirsten HWJ Ten Tusscher, and Alexander V Panfilov.  Effects of heterogeneous diffuse fibrosis on arrhythmia dynamics and mechanism.  Scientific Reports, 6(1):20835, 2016.
  * [23] S Sridhar, Nele Vandersickel, and Alexander V Panfilov.  Effect of myocyte-fibroblast coupling on the onset of pathological dynamics in a model of ventricular tissue.  Scientific Reports, 7(1):40985, 2017.
  * [24] AMGL De Maziere, ACG Van Ginneken, R Wilders, HJ Jongsma, and LN Bouman.  Spatial and functional relationship between myocytes and fibroblasts in the rabbit sinoatrial node.  Journal of Molecular and Cellular Cardiology, 24(6):567–578, 1992.
  * [25] Eva A Rog-Zielinska, Russell A Norris, Peter Kohl, and Roger Markwald.  The living scar–cardiac fibroblasts and the injured heart.  Trends in Molecular Medicine, 22(2):99–114, 2016.
  * [26] Nicola L Walker, Francis L Burton, Sarah Kettlewell, Godfrey L Smith, and Stuart M Cobbe.  Mapping of epicardial activation in a rabbit model of chronic myocardial infarction: response to atrial, endocardial and epicardial pacing.  Journal of Cardiovascular Electrophysiology, 18(8):862–868, 2007.
  * [27] Peter Kohl, Patrizia Camelliti, Francis L Burton, and Godfrey L Smith.  Electrical coupling of fibroblasts and myocytes: relevance for cardiac propagation.  Journal of Electrocardiology, 38(4):45–50, 2005.
  * [28] Kirsten HWJ Ten Tusscher and Alexander V Panfilov.  Influence of diffuse fibrosis on wave propagation in human ventricular tissue.  Europace, 9(suppl_6):vi38–vi45, 2007.
  * [29] Richard H Clayton.  Dispersion of recovery and vulnerability to re-entry in a model of human atrial tissue with simulated diffuse and focal patterns of fibrosis.  Frontiers in Physiology, 9:1052, 2018.
  * [30] David Jakes, Kevin Burrage, Christopher C Drovandi, Pamela Burrage, Alfonso Bueno-Orovio, Rodrigo Weber dos Santos, Blanca Rodriguez, and Brodie AJ Lawson.  Perlin noise generation of physiologically realistic patterns of fibrosis.  BioRxiv, page 668848, 2019.
  * [31] Ron Milo, Shai Shen-Orr, Shalev Itzkovitz, Nadav Kashtan, Dmitri Chklovskii, and Uri Alon.  Network motifs: simple building blocks of complex networks.  Science, 298(5594):824–827, 2002.
  * [32] Sen Song, Per Jesper Sjöström, Markus Reigl, Sacha Nelson, and Dmitri B Chklovskii.  Highly nonrandom features of synaptic connectivity in local cortical circuits.  PLoS Biology, 3(3):e68, 2005.
  * [33] Uri Alon.  Network motifs: theory and experimental approaches.  Nature Reviews Genetics, 8(6):450–461, 2007.
  * [34] Kirsten HWJ ten Tusscher, Denis Noble, Peter-John Noble, and Alexander V Panfilov.  A model for human ventricular tissue.  American Journal of Physiology-Heart and Circulatory Physiology, 286(4):H1573–H1589, 2004.
  * [35] Kirsten HWJ Ten Tusscher and Alexander V Panfilov.  Alternans and spiral breakup in a human ventricular tissue model.  American Journal of Physiology-Heart and Circulatory Physiology, 291(3):H1088–H1100, 2006.
  * [36] K Andrew MacCannell, Hojjat Bazzazi, Lisa Chilton, Yoshiyuki Shibukawa, Robert B Clark, and Wayne R Giles.  A mathematical model of electrotonic interactions between ventricular myocytes and fibroblasts.  Biophysical Journal, 92(11):4121–4132, 2007.
  * [37] L Chilton, S Ohya, D Freed, EDVSY George, V Drobic, Y Shibukawa, KA Maccannell, Y Imaizumi, RB Clark, IMC Dixon, et al.  K+ currents regulate the resting membrane potential, proliferation, and contractile responses in ventricular fibroblasts and myofibroblasts.  American Journal of Physiology-Heart and Circulatory Physiology, 288(6):H2931–H2939, 2005.
  * [38] Alexander Kursanov, Nathalie A Balakina-Vikulova, Olga Solovyova, Alexander Panfilov, and Leonid B Katsnelson.  In silico analysis of the contribution of cardiomyocyte-fibroblast electromechanical interaction to the arrhythmia.  Frontiers in Physiology, 14:390, 2023.
  * [39] Zhilin Qu and Alan Garfinkel.  An advanced algorithm for solving partial differential equation in cardiac conduction.  IEEE transactions on biomedical engineering, 46(9):1166–1168, 1999.
  * [40] Ana Simon-Chica, Eike M Wülfers, and Peter Kohl.  Nonmyocytes as electrophysiological contributors to cardiac excitation and conduction.  American Journal of Physiology-Heart and Circulatory Physiology, 325(3):H475–H491, 2023.
  * [41] Sergio Alonso and Markus Bär.  Reentry near the percolation threshold in a heterogeneous discrete model for cardiac tissue.  Physical Review Letters, 110(15):158101, 2013.
  * [42] Farhad Pashakhanloo and Alexander V Panfilov.  Minimal functional clusters predict the probability of reentry in cardiac fibrotic tissue.  Physical Review Letters, 127(9):098101, 2021.
  * [43] Gabriel Balaban, Caroline Mendonça Costa, Bradley Porter, Brian Halliday, Christopher A Rinaldi, Sanjay Prasad, Gernot Plank, Tevfik F Ismail, and Martin J Bishop.  3d electrophysiological modeling of interstitial fibrosis networks and their role in ventricular arrhythmias in non-ischemic cardiomyopathy.  IEEE Transactions on Biomedical Engineering, 67(11):3125–3133, 2020.
  * [44] Sohail Zahid, Hubert Cochet, Patrick M Boyle, Erica L Schwarz, Kaitlyn N Whyte, Edward J Vigmond, Rémi Dubois, Meleze Hocini, Michel Haïssaguerre, Pierre Jaïs, et al.  Patient-derived models link re-entrant driver localization in atrial fibrillation to fibrosis spatial pattern.  Cardiovascular Research, 110(3):443–454, 2016.
  * [45] Tanmay A Gokhale, Huda Asfour, Shravan Verma, Nenad Bursac, and Craig S Henriquez.  Microheterogeneity-induced conduction slowing and wavefront collisions govern macroscopic conduction behavior: A computational and experimental study.  PLoS Computational Biology, 14(7):e1006276, 2018.
  * [46] Megan E Farquhar, Kevin Burrage, Rodrigo Weber Dos Santos, Alfonso Bueno-Orovio, and Brodie AJ Lawson.  Graph-based homogenisation for modelling cardiac fibrosis.  Journal of Computational Physics, 459:111126, 2022.
  * [47] Brodie AJ Lawson, Rodrigo Weber Dos Santos, Ian W Turner, Alfonso Bueno-Orovio, Pamela Burrage, and Kevin Burrage.  Homogenisation for the monodomain model in the presence of microscopic fibrotic structures.  Communications in Nonlinear Science and Numerical Simulation, 116:106794, 2023.
  * [48] Peter Mortensen, Hao Gao, Godfrey Smith, and Radostin D Simitev.  Addendum: Action potential propagation and block in a model of atrial tissue with myocyte–fibroblast coupling.  Mathematical Medicine and Biology: A Journal of the IMA, 38(3):292–298, 2021.

## Acknowledgements

SS and RHC would like to acknowledge EPSRC EP/T017899/1 The SofTMech
Statistical Emulation and Translation Hub for funding SS. SS and RHC would
like to thank Prof Godfrey Smith and Prof Radostin Simitev for useful
suggestions and comments. SS would also like to thank Prof Sitabhra Sinha for
useful discussions.

## Supplementary Information

Figure 10: $2$-parameter plots for all the $3$ motifs describing the effect of
the coupling strengths $G_{Loc},G_{Long}$ on the APD for both the myocytes
with Shallow parameters. Figure 11: Action potential profiles for both
myocytes in $Motif-1$ while pacing the cells at $T=600$ ms at
$G_{Loc}=G_{Long}=0.5$ nS (a), $2.0$ nS (b) and $4.0$ nS (c) respectively. The
blue and red traces correspond to Shallow and Steep parameter sets. Figure 12:
Action potential profiles for both myocytes in $Motif-3$ while pacing the
cells at $T=600$ ms at $G_{Loc}=G_{Long}=0.5$ nS (a), $2.0$ nS (b) and $4.0$
nS (c) respectively. The blue and red traces correpsond to Shallow and Steep
parameter sets. Figure 13: Effect of delay in stimulation of $myocyte-2$ on
the difference in the APD of the two myocytes for the different motifs with
Shallow parameters. The values of the delay are $\tau_{D}=0$ ms (a,d,g), $=10$
ms (b, e, h)and $=25$ ms (c, f, i). Figure 14: The time series of trans-
membrane potential for both myocytes in $Motif-3$, with $myocyte-2$ stimulated
$25$ ms after $myocyte-1$. Figure 15: $2$-parameter conductance map for pacing
period $T=600$ ms describing the different dynamical regimes for the Steep
parameter set, viz., $NR$, $IR$ and $1:1$, characterised as the fraction of
stimuli that elicit an action potential in the Response Cell for $Motif-1$
(a,b) and $Motif-2$ (c-f). Figure 16: $2$-parameter conductance map for pacing
period $T=600$ ms describing the different dynamical regimes for the Shallow
parameter set, viz., $NR$, $IR$ and $1:1$, characterised as the fraction of
stimuli that elicit an action potential in the Response Cell for $Motif-1$
(a,b) and $Motif-2$ (c-f). Figure 17: $2$-parameter conductance map describing
the different dynamical regimes for the Shallow parameter, viz., $NR$, $IR$
and $1:1$, characterised as the fraction of stimuli that elicit an action
potential in the Response Cell for $Motif-2$ at pacing period $T=300$ ms.
Figure 18: $2$-parameter conductance map describing the different dynamical
regimes for the Steep parameter set, viz., $NR$, $IR$ and $1:1$, characterised
as the fraction of stimuli that elicit an action potential in the Response
Cell for $Motif-2$ at pacing period $T=300$ ms. Figure 19: $2$-parameter
conductance map describing the different dynamical regimes for the Shallow
parameter, viz., $NR$, $IR$ and $1:1$, characterised as the fraction of
stimuli that elicit an action potential in the Response Cell for $Motif-1$
(a,b) and $Motif-3$ (c,d) at pacing period $T=300$ ms. Figure 20:
$2$-parameter conductance map describing the different dynamical regimes for
the Steep parameter set, viz., $NR$, $IR$ and $1:1$, characterised as the
fraction of stimuli that elicit an action potential in the Response Cell for
$Motif-1$ (a,b) and $Motif-3$ (c,d) at pacing period $T=300$ ms. Figure 21:
Effect of coupling motifs on APD for Steep parameters with $V_{FR}=-49$ mV.
Figure 22: Identifying regimes for Steep parameters at $T=600$ ms that
initiate action potentials in non-stimulated cells for $V_{FR}=-49$ mV. Figure
23: Effect of coupling motifs with Shallow cell parameters on APD for
$V_{FR}=-49$ mV. Figure 24: Identifying conductance parameters that initiate
action potentials in non-stimulated myocyte with Shallow restitution for
$V_{FR}=-49$ mV.


2015

Tang

Non-monotone Two-Stage Submodular Maximization

Data Summarization beyond Monotonicity: Non-monotone Two-Stage Submodular
Maximization

Shaojie Tang Naveen Jindal School of Management, The University of Texas at
Dallas

The objective of a two-stage submodular maximization problem is to reduce the
ground set using provided training functions that are submodular, with the aim
of ensuring that optimizing new objective functions over the reduced ground
set yields results comparable to those obtained over the original ground set.
This problem has applications in various domains including data summarization.
Existing studies often assume the monotonicity of the objective function,
whereas our work pioneers the extension of this research to accommodate non-
monotone submodular functions. We have introduced the first constant-factor
approximation algorithms for this more general case.

##  1 Introduction

In this paper, we are motivated by the application of data summarization (Wei
et al. 2013, Mirzasoleiman et al. 2016, Wei et al. 2015, Lin and Bilmes 2011)
and tackle the two-stage submodular optimization problem. In these
applications, we are often faced with multiple user-specific submodular
functions, which are used to evaluate the value of a set of items. A typical
objective is to select a set of $k$ items to maximize each submodular function
(Krause and Golovin 2014). While maximizing a single submodular function has
been widely explored in the literature, the feasibility of existing solutions
diminishes when confronted with a substantial number of submodular functions
and items. Consequently, our objective is to reduce the size of the ground set
in a manner that minimizes the loss when optimizing a new submodular function
over the reduced ground set, as compared to the original ground set.

The problem at hand can be framed as a two-stage submodular maximization
problem, as initially introduced in (Balkanski et al. 2016). While the
majority of prior studies in this domain presume that each submodular function
exhibits monotone non-decreasing behavior, real-world scenarios often involve
objective functions that are non-monotone. These instances include feature
selection (Das and Kempe 2008), profit maximization (Tang and Yuan 2021),
maximum cut (Gotovos et al. 2015), and data summarization (Mirzasoleiman et
al. 2016). A significant contribution presented in our work is the development
of the first constant-factor approximation algorithm for the non-monotone two-
stage submodular maximization problem, with an approximation ratio of $1/2e$.
Remarkably, when the objective function is monotone, our algorithm achieves an
improved approximation ratio of $(1-1/e^{2})/2$, thereby recovering the result
presented in (Stan et al. 2017).

###  1.1 Related Work

The problem of non-monotone submodular maximization has garnered substantial
attention in the literature (Gharan and Vondrák 2011, Buchbinder et al. 2014).
The current state-of-the-art solution for this problem, especially when
accounting for a cardinality constraint, is a $0.385$-approximation algorithm
(Buchbinder and Feldman 2019). However, it is noteworthy that even though each
individual objective function considered in our problem exhibits
submodularity, the overall objective function is not submodular in general. As
a result, the existing findings on non-monotone submodular maximization do not
directly apply to our specific setting.

The most closely related work to our research is the study by (Balkanski et
al. 2016, Mitrovic et al. 2018) and (Stan et al. 2017). They have developed
constant factor approximation algorithms, primarily tailored for the monotone
case. Our work builds upon and extends their results to address the more
general and challenging non-monotone scenario. To achieve this goal, we have
integrated the “local-search” approach (Stan et al. 2017) with “sampling”
technique (Tang 2021) in a non-trivial way, resulting in the creation of a
novel sampling-based algorithm. Furthermore, we have incorporated a trimming
phase into our algorithm, enabling us to attain the first constant-factor
approximation ratio for the non-monotone case.

##  2 Problem Formulation

The input of our problem is a set of $n$ items $\Omega$. There is a group of
$m$ non-monotone submodular functions
$f_{1},\cdots,f_{m}:2^{\Omega}\rightarrow\mathbb{R}_{\geq 0}$. Let
$\Delta_{i}(x,A)=f_{i}(\\{x\\}\cup A)-f_{i}(A)$ denote the marginal gain of
adding $x$ to the set $A$ when considering the function $f_{i}$. Here we say
$f_{i}$ is submodular if and only if
$\Delta_{i}(x,A)\geq\Delta_{i}(x,A^{\prime})$ for any two sets $A$ and
$A^{\prime}$ such that $A\subseteq A^{\prime}\subseteq\Omega$, and any item
$x\in\Omega$ such that $x\notin A^{\prime}$.

Our objective is to compute a reduced ground set $S$ of size $l$, where $l\ll
n$, such that it yields good performance across all $m$ functions when the
choice is limited to items in $S$. Formally, let

$\displaystyle F(S)=\sum_{i\in[m]}\max_{A\subseteq S:|A|\leq k}f_{i}(A)$ (1)

where $k$ is the size constraint of a feasible solution. Our goal is to find
an optimal solution $O\subseteq\Omega$ that maximizes $F$, i.e.,

$\displaystyle O\in\argmax_{S\subseteq\Omega:|S|\leq l}F(S).$ (2)

It is worth mentioning that the objective function $F(\cdot)$ is typically
non-submodular, as observed in (Balkanski et al. 2016). Consequently,
classical algorithms designed for submodular optimization may not provide any
approximation guarantees.

##  3 Algorithm Design and Analysis

Before presenting our algorithm, we need some additional notations. For each
$i\in[m]$, we define the gain associated with removing an item $y$ and
replacing it with $x$ as $\nabla_{i}(x,y,A)=f_{i}(\\{x\\}\cup
A\setminus\\{y\\})-f_{i}(A)$. Then for each $i\in[m]$, we define the largest
possible gain brought by $x$, through local-search, with respect to an
existing set $A$ as $\nabla_{i}(x,A)$. Here the local-search can be realized
either by directly adding $x$ to $A$ (while maintaining the cardinality
constraint) or by substituting it with an item from $T_{i}$. Formally,

$\displaystyle\nabla_{i}(x,A)=\begin{cases}0&\mbox{if }x\in A\\\
\max\\{0,\max_{y\in A}\nabla_{i}(x,y,A),\Delta_{i}(x,A)\\}&\mbox{if }x\notin
A\mbox{ and }|A|<k\\\ \max\\{0,\max_{y\in A}\nabla_{i}(x,y,A)\\}&\mbox{if
}x\notin A\mbox{ and }|A|=k\end{cases}$ (3)

Let $\textsf{Rep}_{i}(x,A)$ represent the item in $A$ that, when replaced by
$x$, maximizes the incremental gain while maintaining feasibility. Formally,

$\displaystyle\textsf{Rep}_{i}(x,A)=\begin{cases}\emptyset&\mbox{if
}\nabla_{i}(x,A)=0\\\ \emptyset&\mbox{if }\nabla_{i}(x,A)>0\mbox{ and
}|A|<k\\\ &\mbox{ and }\max_{y\in A}\nabla_{i}(x,y,A)<\Delta_{i}(x,A)\\\
\arg\max_{y\in A}\nabla_{i}(x,y,A)&\mbox{if }\nabla_{i}(x,A)>0\mbox{ and
}|A|<k\\\ &\mbox{ and }\max_{y\in A}\nabla_{i}(x,y,A)\geq\Delta_{i}(x,A)\\\
\arg\max_{y\in A}\nabla_{i}(x,y,A)&\mbox{if }\nabla_{i}(x,A)>0\mbox{ and
}|A|=k\end{cases}$ (4)

Now we are ready to present the design of our algorithm Sampling-Greedy
(Algorithm 1). Initially, we introduce a set $\Phi$ of $l$ dummy items to the
ground set $\Omega$ to construct an extended ground set
$\Omega^{\prime}=\Omega\cup\Phi$. Here for all $i\in[m]$, we define
$\Delta_{i}(x,A)=0$ for any $x\in\Phi$ and $A\subseteq\Omega^{\prime}$. We
introduce these dummy items to prevent negative marginal utility item
selection. Notably, these dummy items can be safely removed from the solution
without altering its utility.

Throughout the process, Sampling-Greedy maintains a solution set denoted as
$S$, along with a set of feasible solutions $T_{i}\subseteq S$ for each
function $f_{i}$ (all of which are initially set to empty). In each iteration,
it first computes the top $l$ items $M$ from the extended ground set
$\Omega^{\prime}$ based on its combined contribution to each $f_{i}$,
indicated by $\sum_{i=1}^{m}\nabla_{i}(x,T_{i})$. That is,

$\displaystyle M=\argmax_{A\subseteq\Omega^{\prime}:|A|=l}\sum_{x\in
A}\sum_{i=1}^{m}\nabla_{i}(x,T_{i}).$ (5)

Then it randomly selects one item, say $x^{*}$, from $M$ and adds $x^{*}$ to
$S$. Sampling-Greedy then verifies if any of the sets $T_{i}$ can be improved.
This can be achieved by either directly adding $x^{*}$ (while adhering to the
cardinality constraint) or substituting it with an item from $T_{i}$. For each
$i\in[m]$, we update $T_{i}$ if and only if $\nabla_{i}(x^{*},A)>0$.

Note that there might exist some $i\in[m]$ and $x\in T_{i}$ such that
$f_{i}(T_{i})-f_{i}(T_{i}\setminus\\{x\\})<0$. In other words, certain subsets
$T_{i}$ could contain items that provide negative marginal utility to the set
$T_{i}$. Consequently, we introduce a “trimming” phase (Algorithm 2) to refine
each $T_{i}$ and ensure that no item contributes negative utility to it. This
can be achieved through an iterative process of evaluating the marginal
utility of each item within $T_{i}$ and subsequently removing any items with
negative marginal utility. By the submodularity of $f_{i}$, we can show that
after this trimming phase, $T_{i}$ does not contain any items whose marginal
utility if negative. It is also easy to verify that the trimming phase does
not decrease the utility of our solution. A formal description of these
properties is presented in the following lemma.

Algorithm 1 Sampling-Greedy

1: $S\leftarrow\emptyset,T_{i}\leftarrow\emptyset(\forall i\in[m])$

2: for $j\in[l]$ do

3: $M=\argmax_{A\subseteq\Omega^{\prime}:|A|=l}\sum_{x\in
A}\sum_{i=1}^{m}\nabla_{i}(x,T_{i}).$

4: randomly pick one item $x^{*}$ from $M$, $S\leftarrow S\cup\\{x^{*}\\}$

5: for $i\in[m]$ do

6: if $\nabla_{i}(x^{*},A)>0$ then

7: $T_{i}\leftarrow
T_{i}\setminus\textsf{Rep}_{i}(x^{*},T_{i})\cup\\{x^{*}\\}$

8: $T_{i}\leftarrow\textsf{Trim}(T_{i},f_{i})$

9: return $S,T_{1},T_{2},\cdots,T_{m}$

Algorithm 2 Trim($B$, $f_{i}$)

1: $A\leftarrow B$

2: for $x\in A$ do

3: if $f_{i}(A)-f_{i}(A\setminus\\{x\\})<0$ then

4: $A\leftarrow A\setminus\\{x\\}$

5: return $A$

###### Lemma 3.1

Consider any set of items $B\subseteq\Omega$ and a function $f_{i}$. Assume
$A$ is returned from Trim($B$, $f_{i}$), we have $f_{i}(A)\geq f_{i}(B)$ and
for all $x\in A$, we have $f_{i}(A)-f_{i}(A\setminus\\{x\\})\geq 0$.

_Proof:_ The proof that $f_{i}(A)\geq f_{i}(B)$ is straightforward, as it
follows from the fact that the trimming phase only eliminates items with a
negative marginal contribution. We next prove that for all $x\in A$, we have
$f_{i}(A)-f_{i}(A\setminus\\{x\\})\geq 0$. We prove this through
contradiction. Suppose there exists an item $y\in A$ such that
$f_{i}(A)-f_{i}(A\setminus\\{y\\})<0$. Let’s denote the solution before
considering the inclusion of $y$ as $A^{\prime}$. In this case, it must hold
that $f_{i}(A^{\prime})-f_{i}(A^{\prime}\setminus\\{y\\})\geq 0$, as
otherwise, the trimming phase would eliminate $y$ from the solution.
Furthermore, it is straightforward to confirm that $A\subseteq A^{\prime}$. As
a consequence, based on the assumption that $f_{i}$ is a submodular function,
we have $f_{i}(A)-f_{i}(A\setminus\\{y\\})\geq
f_{i}(A^{\prime})-f_{i}(A^{\prime}\setminus\\{y\\})$. This, together with
$f_{i}(A^{\prime})-f_{i}(A^{\prime}\setminus\\{y\\})\geq 0$, implies that
$f_{i}(A)-f_{i}(A\setminus\\{y\\})\geq
f_{i}(A^{\prime})-f_{i}(A^{\prime}\setminus\\{y\\})\geq 0$. This contradicts
to the assumption that $f_{i}(A)-f_{i}(A\setminus\\{y\\})<0$. $\Box$

###  3.1 Performance analysis

First, it is easy to verify that Sampling-Greedy requires $O(l(mkl+mn))$
function evaluations. This is because Sampling-Greedy comprises $l$
iterations, where each iteration involves $mkl$ function evaluations in Line 3
of Algorithm 1, along with an additional $mn$ function evaluations in
Algorithm 2. In the following theorem, we show that the expected utility of
our solution is at least a constant factor approximation of the optimal
solution.

###### Theorem 3.2

Sampling-Greedy returns a random set $S$ of size at most $l$ such that

$\displaystyle\mathbb{E}_{S}[F(S)]\geq\frac{1}{2e}F(O)$ (6)

where $O$ represents the optimal solution.

The rest of this section is devoted to proving this theorem. The basic idea
behind the proof is to establish a lower bound on the expected marginal
utility achieved by adding $x^{*}$ to set $S$ after each iteration. We
demonstrate that this utility increment is substantial enough to guarantee a
$1/2e$ approximation ratio. Consider an arbitrary round $t\in[l]$ of Sampling-
Greedy, let $S$ and $T_{1},\cdots,T_{m}$ denote the solution obtained at the
end of round $t$. By the design of Sampling-Greedy, we randomly pick an item
$x^{*}$ from $M$ and add it to $S$, hence, by the definition of $M$, the
expected marginal utility of adding $x^{*}$ to $S$ before the “trimming phase”
is

$\displaystyle\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]=\frac{1}{l}\max_{A\subseteq\Omega^{\prime}:|A|=l}\sum_{x\in
A}\sum_{i=1}^{m}\nabla_{i}(x,T_{i}).$ (7)

Recall that the trimming phase does not decrease utility. Therefore, the
ultimate expected utility increment after each iteration is at least
$\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]$. Moreover, because
$F$ is a monotone function, it is safe to assume that the size of the optimal
solution is $l$, i.e, $|O|=l$. We next provide a lower bound on
$\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]$.

Observe that

$\displaystyle\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]=\frac{1}{l}\max_{A\subseteq\Omega^{\prime}:|A|=l}\sum_{x\in
A}\sum_{i=1}^{m}\nabla_{i}(x,T_{i})\geq\frac{1}{|O|}\sum_{x\in
O}\sum_{i=1}^{m}\nabla_{i}(x,T_{i})=\frac{1}{l}\sum_{x\in
O}\sum_{i=1}^{m}\nabla_{i}(x,T_{i})$ (8)

Let $O_{i}\subseteq O$ represent a subset with a maximum size of $k$ items,
chosen to maximize $f_{i}$, i.e., $O_{i}=\argmax_{A\subseteq O:|A|\leq
k}f_{i}(A)$. Inequality (8) implies that

$\displaystyle\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]\geq\frac{1}{l}\sum_{x\in
O}\sum_{i=1}^{m}\nabla_{i}(x,T_{i})\geq\frac{1}{l}\sum_{i=1}^{m}\sum_{x\in
O_{i}}\nabla_{i}(x,T_{i}).$ (9)

For the sake of simplifying our analysis, let’s assume that the size of
$O_{i}$ is fixed at $k$, i.e, $|O_{i}|=k$. If $O_{i}$ does not initially
contain exactly $k$ items, we can add additional dummy items to bring its size
up to $k$. Please note that if $x$ is a dummy item, then
$\nabla_{i}(x,T_{i})=0$. This stems from the observation that $T_{i}$ does not
contain any items that would yield negative marginal utility, a consequence of
the trimming phase (Lemma 3.1). Hence, adding dummy items to $O_{i}$ does not
affect inequality (9).

It is easy to verify that there is a mapping $\pi$ between $O_{i}$ and $T_{i}$
such that every item of $O_{i}\cap T_{i}$ is mapped to itself, and every item
of $O_{i}\setminus T_{i}$ is mapped to either the empty set or an item in
$T_{i}\setminus O_{i}$. We next give a lower bound of $\nabla_{i}(x,T_{i})$.

###### Lemma 3.3

For all $i\in[m]$ and $x\in O_{i}$, we have

$\displaystyle\nabla_{i}(x,T_{i})\geq\Delta_{i}(x,T_{i})-\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\}).$
(10)

_Proof:_ We prove this lemma in three cases. We first consider the case when
$x\notin T_{i}$ and $\pi(x)\neq\emptyset$. In this case, the following chain
proves this lemma.

$\displaystyle\nabla_{i}(x,T_{i})$ $\displaystyle\geq f_{i}(\\{x\\}\cup
T_{i}\setminus\\{\pi(x)\\})-f_{i}(T_{i})$ (13)
$\displaystyle=\Delta_{i}(x,T_{i})-\Delta_{i}(\pi(x),T_{i}\cup\\{x\\}\setminus\\{\pi(x)\\})$
$\displaystyle\geq\Delta_{i}(x,T_{i})-\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\})$

where the first inequality is by the definition of $\nabla_{i}(x,T_{i})$ and
the second inequality is by the assumption that $f_{i}$ is a submodular
function.

We next consider the case when $x\notin T_{i}$ and $\pi(x)=\emptyset$. In this
case, because $\pi(x)=\emptyset$, i.e., $x$ is not mapped to any item from
$T_{i}$, we have $|T_{i}|<k$. Hence,

$\displaystyle\nabla_{i}(x,T_{i})=\max\\{0,\max_{y\in
T_{i}}\nabla_{i}(x,y,T_{i}),\Delta_{i}(x,T_{i})\\}\geq\Delta_{i}(x,T_{i}).$
(14)

Moreover, $\pi(x)=\emptyset$ implies that

$\displaystyle\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\})=0.$ (15)

It follows that

$\displaystyle\nabla_{i}(x,T_{i})\geq\Delta_{i}(x,T_{i})-0=\Delta_{i}(x,T_{i})-\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\}),$
(16)

where the inequality is by inequality (14) and the equality is by equality
(15).

At last, we consider the case when $x\in T_{i}$. In this case, we have
$\Delta_{i}(x,T_{i})=0$, and
$\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\})\geq 0$, a consequence of the
trimming phase (Lemma 3.1). Hence,
$\Delta_{i}(x,T_{i})-\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\})\leq 0$. It
follows that

$\displaystyle\nabla_{i}(x,T_{i})\geq
0\geq\Delta_{i}(x,T_{i})-\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\}).$ (17)

$\Box$

Inequality (9) and Lemma 3.3 imply that

$\displaystyle\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]$
$\displaystyle\geq\frac{1}{l}\sum_{i=1}^{m}\sum_{x\in
O_{i}}\nabla_{i}(x,T_{i})$ (19)
$\displaystyle\geq\frac{1}{l}\sum_{i=1}^{m}\sum_{x\in
O_{i}}(\Delta_{i}(x,T_{i})-\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\})).$

Because $f_{i}$ is submodular, we have

$\displaystyle\sum_{x\in O_{i}}\Delta_{i}(x,T_{i})\geq f_{i}(O_{i}\cup
T_{i})-f_{i}(T_{i}).$ (20)

Moreover, no two items from $O_{i}$ are mapped to the same item from $T_{i}$,
we have

$\displaystyle\sum_{x\in
O_{i}}\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\})=\sum_{y\in
T_{i}}\Delta_{i}(y,T_{i}\setminus\\{y\\})\leq f_{i}(T_{i})$ (21)

where the inequality is by the assumption that $f_{i}$ is submodular.

Inequalities (19), (20) and (21) together imply that

$\displaystyle\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]$
$\displaystyle\geq\frac{1}{l}\sum_{i=1}^{m}\sum_{x\in
O_{i}}(\Delta_{i}(x,T_{i})-\Delta_{i}(\pi(x),T_{i}\setminus\\{\pi(x)\\}))$
(24) $\displaystyle\geq\frac{1}{l}\sum_{i=1}^{m}(f_{i}(O_{i}\cup
T_{i})-f_{i}(T_{i})-f_{i}(T_{i}))$
$\displaystyle=\frac{1}{l}\sum_{i=1}^{m}(f_{i}(O_{i}\cup
T_{i})-2f_{i}(T_{i})).$

Taking the expectation over $T_{1},\cdots,T_{m}$ for both the left and right
hand sides of (24), we have

$\displaystyle\mathbb{E}_{T_{1},\cdots,T_{m}}\big{[}\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]\big{]}$
$\displaystyle\geq\mathbb{E}_{T_{1},\cdots,T_{m}}[\frac{1}{l}\sum_{i=1}^{m}(f_{i}(O_{i}\cup
T_{i})-2f_{i}(T_{i}))]$ (29)
$\displaystyle=\mathbb{E}_{T_{1},\cdots,T_{m}}[\frac{1}{l}\sum_{i=1}^{m}(f_{i}(O_{i}\cup
T_{i}))]-\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}\frac{2}{l}f_{i}(T_{i}))]$
$\displaystyle=\frac{1}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}(f_{i}(O_{i}\cup
T_{i}))]-\frac{2}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}f_{i}(T_{i}))]$
$\displaystyle\geq\frac{1}{l}(1-\frac{1}{l})^{t}\sum_{i=1}^{m}f_{i}(O_{i})-\frac{2}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[f_{i}(T_{i}))]$
$\displaystyle=\frac{1}{l}(1-\frac{1}{l})^{t}F(O)-\frac{2}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[f_{i}(T_{i}))].$

The second inequality is by the observation that
$\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}(f_{i}(O_{i}\cup
T_{i}))]\geq(1-\frac{1}{l})^{t}\sum_{i=1}^{m}f_{i}(O_{i})$. To prove this
inequality, recall that in each round, Sampling-Greedy randomly picks an item
from $M$ to be included in $S$. Hence, right before entering round $t$ of
Sampling-Greedy, each item $x\in\Omega^{\prime}$ has a probability of at most
$p=1-(1-\frac{1}{l})^{t}$ of being included in $S$ and consequently in $T_{i}$
for all $i\in[m]$. By Lemma 2.2 of (Buchbinder et al. 2014), we have
$\mathbb{E}_{T_{i}}[f_{i}(O_{i}\cup
T_{i})]\geq(1-p)f_{i}(O_{i})=(1-\frac{1}{l})^{t}f_{i}(O_{i})$ for all
$i\in[m]$. It follows that
$\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}(f_{i}(O_{i}\cup
T_{i}))]\geq(1-\frac{1}{l})^{t}\sum_{i=1}^{m}f_{i}(O_{i})$.

Let $X_{t}$ denote the value of
$\mathbb{E}_{T_{1},\cdots,T_{m}}\big{[}\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]\big{]}$
at the end of round $t$. Inequality (29) implies that

$\displaystyle
X_{t+1}-X_{t}\geq\frac{1}{l}(1-\frac{1}{l})^{t+1}F(O)-\frac{2}{l}X_{t}$ (30)
$\displaystyle\Rightarrow
2(X_{t+1}-X_{t})\geq\frac{1}{l}(1-\frac{1}{l})^{t+1}F(O)-\frac{2}{l}X_{t}$
(31) $\displaystyle\Rightarrow
2X_{t+1}-2X_{t}\geq\frac{1}{l}(1-\frac{1}{l})^{t+1}F(O)-\frac{2}{l}X_{t}$ (32)
$\displaystyle\Rightarrow
2X_{t+1}\geq\frac{1}{l}(1-\frac{1}{l})^{t+1}F(O)+(2-\frac{2}{l})X_{t}.$ (33)

Based on the above inequality, we next prove through induction that
$2X_{t}\geq\frac{t}{l}(1-\frac{1}{l})^{t-1}F(O)$. Note that $X_{0}=0$, meaning
that the utility before the start of the algorithm is zero. The induction step
is established in the following manner:

$\displaystyle
2X_{t+1}\geq\frac{1}{l}(1-\frac{1}{l})^{t+1}F(O)+(2-\frac{2}{l})X_{t}$ (34)
$\displaystyle\Rightarrow 2X_{t+1}$
$\displaystyle\geq\frac{1}{l}(1-\frac{1}{l})^{t+1}F(O)+(2-\frac{2}{l})\frac{t}{l}(1-\frac{1}{l})^{t-1}F(O)$
(40)
$\displaystyle=\frac{1}{l}(1-\frac{1}{l})^{t+1}F(O)+2\frac{t}{l}(1-\frac{1}{l})^{t}F(O)$
$\displaystyle=(\frac{1}{l}(1-\frac{1}{l})+2\frac{t}{l})(1-\frac{1}{l})^{t}F(O)$
$\displaystyle\geq 2\frac{t}{l}(1-\frac{1}{l})^{t}F(O)$
$\displaystyle=(\frac{t}{l}+\frac{t}{l})(1-\frac{1}{l})^{t}F(O)$
$\displaystyle\geq\frac{t+1}{l}(1-\frac{1}{l})^{t}F(O).$

It follows that the value of $2X_{l}$ is at least $(1-\frac{1}{l})^{l-1}F(O)$,
which itself is bounded from below by $(1/e)\cdot F(O)$. Here, $X_{l}$
represents the expected utility of our algorithm upon completion. Hence, the
expected utility of our algorithm is at least $X_{l}\geq(1/2e)\cdot F(O)$.

###  3.2 Enhanced results for monotone case

For the case when $f_{i}$ is both monotone and submodular, we will demonstrate
that the approximation ratio of Sampling-Greedy is improved to $(1-1/e^{2})/2$
which recovers the results presented in (Stan et al. 2017). Observe that if
$f_{i}$ is monotone, we have $f_{i}(O_{i}\cup T_{i})\geq f_{i}(O_{i})$. Hence,
inequality (29) implies that

$\displaystyle\mathbb{E}_{T_{1},\cdots,T_{m}}\big{[}\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]\big{]}$
$\displaystyle\geq\frac{1}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}(f_{i}(O_{i}\cup
T_{i}))]-\frac{2}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}f_{i}(T_{i}))]$
(44)
$\displaystyle\geq\frac{1}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}(f_{i}(O_{i}))]-\frac{2}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}f_{i}(T_{i}))]$
$\displaystyle=\frac{1}{l}\sum_{i=1}^{m}f_{i}(O_{i})-\frac{2}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}f_{i}(T_{i}))]$
$\displaystyle=\frac{1}{l}F(O)-\frac{2}{l}\mathbb{E}_{T_{1},\cdots,T_{m}}[\sum_{i=1}^{m}f_{i}(T_{i}))]$

where the first equality is because $O_{i}$ is a fixed set for all $i\in[m]$.
Let $X_{t}$ denote the value of
$\mathbb{E}_{T_{1},\cdots,T_{m}}\big{[}\mathbb{E}_{x^{*}}[\sum_{i=1}^{m}\nabla_{i}(x^{*},T_{i})]\big{]}$
at the end of round $t$. Inequality (44) implies that

$\displaystyle X_{t+1}-X_{t}\geq\frac{1}{l}F(O)-\frac{2}{l}X_{t}.$ (45)

Previous research (Stan et al. 2017) has demonstrated that by inductively
solving the equation above, we can establish that
$X_{l}\geq((1-1/e^{2})/2)\cdot F(O)$.

## References

  * Balkanski et al. (2016) Balkanski, Eric, Baharan Mirzasoleiman, Andreas Krause, Yaron Singer. 2016.  Learning sparse combinatorial representations via two-stage submodular maximization.  International Conference on Machine Learning. PMLR, 2207–2216.
  * Buchbinder and Feldman (2019) Buchbinder, Niv, Moran Feldman. 2019.  Constrained submodular maximization via a nonsymmetric technique.  Mathematics of Operations Research 44 988–1005.
  * Buchbinder et al. (2014) Buchbinder, Niv, Moran Feldman, Joseph Naor, Roy Schwartz. 2014.  Submodular maximization with cardinality constraints.  Proceedings of the twenty-fifth annual ACM-SIAM symposium on Discrete algorithms. SIAM, 1433–1452.
  * Das and Kempe (2008) Das, Abhimanyu, David Kempe. 2008.  Algorithms for subset selection in linear regression.  Proceedings of the fortieth annual ACM symposium on Theory of computing. 45–54.
  * Gharan and Vondrák (2011) Gharan, Shayan Oveis, Jan Vondrák. 2011.  Submodular maximization by simulated annealing.  Proceedings of the twenty-second annual ACM-SIAM symposium on Discrete Algorithms. SIAM, 1098–1116.
  * Gotovos et al. (2015) Gotovos, Alkis, Amin Karbasi, Andreas Krause. 2015.  Non-monotone adaptive submodular maximization.  Twenty-Fourth International Joint Conference on Artificial Intelligence.
  * Krause and Golovin (2014) Krause, Andreas, Daniel Golovin. 2014.  Submodular function maximization.  Tractability 3 3\.
  * Lin and Bilmes (2011) Lin, Hui, Jeff Bilmes. 2011.  A class of submodular functions for document summarization.  Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies. 510–520.
  * Mirzasoleiman et al. (2016) Mirzasoleiman, Baharan, Ashwinkumar Badanidiyuru, Amin Karbasi. 2016.  Fast constrained submodular maximization: Personalized data summarization.  ICML. 1358–1367.
  * Mitrovic et al. (2018) Mitrovic, Marko, Ehsan Kazemi, Morteza Zadimoghaddam, Amin Karbasi. 2018.  Data summarization at scale: A two-stage submodular approach.  International Conference on Machine Learning. PMLR, 3596–3605.
  * Stan et al. (2017) Stan, Serban, Morteza Zadimoghaddam, Andreas Krause, Amin Karbasi. 2017.  Probabilistic submodular maximization in sub-linear time.  International Conference on Machine Learning. PMLR, 3241–3250.
  * Tang (2021) Tang, Shaojie. 2021.  Beyond pointwise submodularity: Non-monotone adaptive submodular maximization in linear time.  Theoretical Computer Science 850 249–261.
  * Tang and Yuan (2021) Tang, Shaojie, Jing Yuan. 2021.  Adaptive regularized submodular maximization.  32nd International Symposium on Algorithms and Computation (ISAAC 2021). Schloss Dagstuhl-Leibniz-Zentrum für Informatik.
  * Wei et al. (2015) Wei, Kai, Rishabh Iyer, Jeff Bilmes. 2015.  Submodularity in data subset selection and active learning.  International conference on machine learning. PMLR, 1954–1963.
  * Wei et al. (2013) Wei, Kai, Yuzong Liu, Katrin Kirchhoff, Jeff Bilmes. 2013.  Using document summarization techniques for speech data subset selection.  Proceedings of the 2013 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies. 721–726.