# Full-frequency dynamic convolution: a physical frequency-dependent
convolution for sound event detection

###### Abstract

Recently, 2D convolution has been found unqualified in sound event detection
(SED). It enforces translation equivariance on sound events along frequency
axis, which is not a shift-invariant dimension. To address this issue, dynamic
convolution is used to model the frequency dependency of sound events. In this
paper, we proposed the first full-dynamic method named _full-frequency dynamic
convolution_ (FFDConv). FFDConv generates frequency kernels for every
frequency band, which is designed directly in the structure for frequency-
dependent modeling. It physically furnished 2D convolution with the capability
of frequency-dependent modeling. FFDConv outperforms not only the baseline by
6.6% in DESED real validation dataset in terms of PSDS1, but outperforms the
other full-dynamic methods. In addition, by visualizing features of sound
events, we observed that FFDConv could effectively extract coherent features
in specific frequency bands, consistent with the vocal continuity of sound
events. This proves that FFDConv has great frequency-dependent perception
ability. Code is available at FFDConv.

Index Terms—  sound event detection, full-frequency dynamic convolution,
frequency-dependent modeling, independent representation spaces, vocal
continuity

##  1 introduction

Sound event detection (SED) is one of the subtasks of computational auditory
scene analysis (CASA) [1], which helps machines understand the content of an
audio scene. Similar to visual object detection [2] and segmentation [3], SED
aims to detect sound events and corresponding timestamps (onset and offset),
considered as a prior task of automatic speech recognition (ASR) and speaker
verification. It has wide applications in information retrieval [4], smart
homes [5], and smart cities [6].

Fig. 1: Illustration of frequency-dependent modeling. Top models time-
frequency patterns in the same space with a shared kernel. Bottom models them
in serval spaces with frequency-adaptive kernels, in which time-frequency
patterns specific to sound events can be considered.

SED has achieved great success with the help of deep learning (DL). The
general paradigm is that acoustic spectral features are passed through a deep
neural network and then transformed into discriminative acoustic
representations to distinguish different sound events. Designing an effective
feature extractor has become a hot topic in SED, which has been adopting
methods qualified in other domains in the past few years. Convolutional neural
networks (CNN) from the field of computer vision, such as SENet [7], SKNet
[8], and CBAM [9] have been migrated to SED in the spirit of acoustic
spectrogram being similar to two-dimensional image data. With the intention
that speech and audio are both sound data, Conformer [10] from the field of
speech recognition has migrated to SED. However, they all failed to show good
performance in SED. Specifically, SENet, SKNet, and CBAM are designed on image
data with a clear 2D spatial concept, while audio data is a sequence.
Conformer is designed on speech data containing only the speech sound event,
meaning time-frequency patterns of speech data are distributed only in a
certain fixed frequency band. However, audio data always contains multiple
sound events, and so has diverse time-frequency patterns of sound events. All
of the above emphasize that DL methods qualified in other domains may not
necessarily be compatible with SED.

Dynamic convolution network [11] was initially proposed for video prediction.
It was designed to generate future frames based on the motion pattern within a
particular video. The parameters of the dynamic convolution kernel are always
adapted to the input. In SED, different sound events are distributed in
different frequency regions, and this frequency dependence is invariant over
time. This has motivated some researchers to investigate whether dynamic
convolution can improve the capability of 2D convolution in modeling the
frequency dependence of sound events. [12] proposed frequency dynamic
convolution (FDConv), which found that the time-frequency spectrogram is not
translation invariant on frequency dimension like image data. FDConv extracts
frequency-adaptive attention weights from input for several pre-initialized
convolution kernels. These kernels are then weightedly combined in the number
dimension to obtain one convolution kernel. Then, the combined kernel is
convoluted with the input in a standard manner. [12, 13], [14] proposed multi-
dimensional frequency dynamic convolution (MFDConv), which extends the
frequency-adaptive dynamic properties of convolutional kernels to more
dimensions of the kernel space, i.e. in-channels, out-channels, and kernel
numbers.

Although FDConv and MFDConv have achieved great performance, they are
essentially the same as basic convolution, which is spatially shared. They
belong to semi-dynamic convolution in the field of dynamic convolution. As
shown in the upper part of Fig. 1, their perception abilities of different
frequency bands are identical. They can only model time-frequency patterns in
one representation space, where sound events are not easily recognized from
each other. Compared with semi-dynamic convolution, full-dynamic convolution
[11, 15, 16, 17, 18] attracts more attention recently, which uses a separate
network branch to predict a specific filter for each pixel. [18] found this
type of dynamic convolution is equivalent to applying attention on unfolded
input features, which enables it more effective when modeling complex
patterns. Sound events’ time-frequency patterns are highly frequency-
dependent, and full-dynamic convolution can model features of spatial pixels
with different filters. Full-dynamic convolution may be optimal in dealing
with recognizing sound events.

In this paper, we propose a novel method named _full-frequency dynamic
convolution_ (FFDConv), which is the first full-dynamic convolution method for
SED. As shown in the lower part of Fig. 1, FFDConv generates frequency-
specific kernels, resulting in distinct representation spaces. This design is
applied directly in the structure for frequency-dependent modeling. In this
way, the 2D convolution is physically furnished with the capability of
frequency-dependent modeling, so that the specific time-frequency patterns can
be acquired for different sound events. In the end, sound events can be easily
recognized from each other in subsequent classification.

The main contributions of this paper are summarized as follows:

  * •

We proposed full-frequency dynamic convolution that can model time-frequency
patterns in independent representation spaces. This method will extract more
discriminative features of sound events, resulting in effective
classification.

  * •

The Proposed method outperforms not only baseline but also pre-existing full
dynamic filters method in other domain.

  * •

By visualizing features of sound events, we found the ability to model
temporally coherent features is essential to the detection of sound events.
And the FFDConv has this ability.

##  2 mthodology

###  2.1 Full-dynamic convolution

A basic 2D convolution can be denoted as $y=\boldsymbol{W}\ast
x+\boldsymbol{b}$, where $x\in{\mathbb{R}^{T\times F\times C_{in}}}$ and
$y\in{\mathbb{R}^{T\times F\times C_{out}}}$ denote the input feature and
output feature; $\boldsymbol{W}\in{\mathbb{R}^{k\times k\times C_{in}\times
C_{out}}}$ and $\boldsymbol{b}\in{\mathbb{R}^{C_{out}}}$ denote the weight and
bias of a basic convolution kernel. In contrast to basis convolution, full-
dynamic convolution [11] leverages separate network branches to generate the
filters for each pixel. Full-dynamic convolution operation can be written as:

$\displaystyle y=\boldsymbol{Concat}$ $\displaystyle{(\boldsymbol{W}_{t,f}\ast
x(t,f))}$ (1) $\displaystyle\boldsymbol{W}_{t,f}=$ $\displaystyle G(x,t,f)$

where $\boldsymbol{W}_{t,f}$ denotes weight for the current pixel; The $G$ is
the filter generating function; $Concat$ here aims to convey that convolution
operation of each pixel is independent. For simplicity, the bias term is
omitted.

Fig. 2: Illustration of full-frequency dynamic convolution. In general, the
factory produces frequency-dependent kernels from acoustic feature, and then
kernels are convoluted with input along the time axis. In the factory, there
are two workshops aiming to produce spatial filters and channel filters,
respectively. And they are integrated in the assembly workshop.

###  2.2 Overall of proposed method

As is commonly understood, different sound events have different frequency
band distributions. For instance, catcall, which is sharp, shrill, and high-
pitched, is often heard in the high-frequency range; running water, which is
low, soft, and soothing, is often heard in the low-frequency range. Based on
this, we explore designing a new convolution for SED, which can capture the
distribution of frequency bands and model time-frequency patterns of sound
events in different frequency representation spaces.

Inspired by full dynamic convolution [18], we designed the full-frequency
dynamic convolution (FFDConv) for SED. Overall, as shown in Fig. 2, FFDConv
employs a separate branch to predict kernels for each frequency band, in which
the content of kernels is based on input feature. In the kernel-generating
branch, there are two sub-branches: the spatial filter-generating branch for
the spatial space of kernels and the channel filter-generating branch for the
channel space of kernels. After spatial and channel filters are obtained, they
are combined and then convoluted with the input feature. Note that similarly,
full-temporal dynamic convolution (FTDConv) predicts kernels for each temporal
frame, and kernels are convoluted with input along the frequency axis.

###  2.3 Full-frequency dynamic convolution

Unlike the previous semi-dynamic convolution, FFDConv is designed directly in
the structure for frequency-dependent modeling. It models the feature along
the frequency axis in different representation spaces. Mathematically, FFDConv
can be written as:

$\displaystyle y=\boldsymbol{Concat}$ $\displaystyle{(\boldsymbol{W}_{f}\ast
x(f),\,dim=f)}$ (2) $\displaystyle\boldsymbol{W}_{f}=$
$\displaystyle\;G_{s}(x,f)\odot G_{c}(x,f)$

where $\boldsymbol{W}_{f}$ is the content-adaptive kernel for the $f^{th}$
frequency band; $x(f)\in\mathbb{R}^{T}$ is the $f^{th}$ frequency band of
input feature; $G_{s}$ and $G_{c}$ are the spatial and channel filter-
generating function; $\odot$ denotes the elemental dot product operator. For
clarity, $Concat$ here aims to convey that $\boldsymbol{W}_{f}$ is convoluted
with input along the time axis.

FFDConv employs a separate branch to generate convolution kernels for each
frequency band, in which there are two sub-branches: spatial filter-generating
branch and channel filter-generating branch. The spatial filter-generating
module is designed to predict the spatial content of dynamic kernels, and the
channel-generating module is designed to predict the channel content of
dynamic kernels. For efficiency, the dynamic filters are decoupled into
spatial and channel ones, following [18].

Spatial filter generating. As illustrated in Fig. 3, we use a standard Conv2D
to compress the time dimension of input and map channel dimension from $C$ to
$K^{2}$, whose kernel weight $W\in\mathbf{R}^{C\times K^{2}\times T\times W}$,
where $W$ is the window size of the kernel in the frequency dimension. It
moves along the frequency axis when convoluted with input. In this way, not
only are the adjacent frequency components considered, but information along
the time axis is aggregated. Then, the spatial filter of FFDConv is obtained,
which assigns $K\times K$ spatial weight to every frequency kernel and is
highly related to the input. Note that full dynamic convolution [18] assigns
$K\times K$ spatial weight to every pixel. Consequently, FFDConv can model
features from different frequency bands of the input in independent
representation spaces.

Considering these representation spaces may be far apart from each other, we
employ an attention module following [12] to limit individual differences
between them so as not to be too large. Finally, the spatial filter is passed
through a Filter-Norm module following [18], avoiding the gradient
vanishing/exploding during training.

Channel filter generating. As illustrated in Fig. 3, the channel filter
generating module is similar to the SE block [7]. It compresses the time and
frequency feature of input by applying an average pooling and maps the channel
dimension from $C$ to $CK^{2}$ by two fully connected (FC) layers. Between two
fully connected layers, the ReLU activation function is applied to introduce
non-linearity. After input is passed through this module, the channel filter
of FFDConv is obtained, which assigns $C$ channel weight to each spatial
location of the frequency kernel. It should be noted that the channel filter
for $F$ frequency kernels is the same. In the end, the channel filter is also
passed through the Filter-Norm [18]. The spatial and channel filters are mixed
by dot product, and the full frequency kernels are obtained. We then use them
to model time-frequency patterns of input features.

###  2.4 FFDConv block

Considering frequency kernels of FFDConv don’t have output channel dimension,
we design an FFDConv block that contains the channel mapping. As illustrated
in Fig. 3, firstly, the channel dimension of input is mapped from $C_{in}$ to
$C_{out}$ after passing through the channel transformation module. Then, based
on the input feature, the spatial and channel filters are obtained by passing
through the spatial and channel filter generating module. Full-frequency
dynamic kernels are obtained by mixing the spatial and channel filters.
Finally, the kernels are convoluted with input along the time axis.

In the actual algorithm, following [18], spatial filters, channel filters, and
input are sent to DDF operation to get the output, which is implemented in
CUDA, alleviating any need to save intermediate multiplied filters during
network training and inference. Note that the DDF op needs $H\times W$ spatial
filters. We repeat the $1\times F$ spatial filters to $T\times F$ so that the
kernel’s weights are the same along the time axis when convoluted with input
in $f^{th}$ frequency band.

Fig. 3: Details of the FFDConv

##  3 experiment

###  3.1 Dataset and experiment setup

All experiments are conducted on the dataset of Task 4 in the DCASE 2022. The
training set consists of three types of data: weakly labeled data (1578
clips), synthetic strongly labeled data (10000 clips), and unlabeled in-domain
data (14412 clips). The real validation set (1168 clips) is used for
evaluation. The input acoustic feature is the log Mel spectrogram extracted
from 10-second-long audio data with a sampling rate of 16 kHz. The feature
configuration is the same as [13], in which the input feature has 626 frames
and 128 mel frequency bands.

The baseline model is the CRNN architecture [19], which consists of 7 layers
of conv blocks and 2 layers of Bi-GRU. Attention pooling module is added at
the last FC layer for joint training of weakly labeled data, and mean teacher
(MT) [20] is applied for consistency training with unlabeled data for semi-
supervised learning. Data augmentations such as MixUp [21], time masking [22],
frame-shift, and FilterAugment [23] are used. The data augmentation parameters
are identical to [12].

Poly-phonic sound event detection scores (PSDS), collar-based F1 score
(EB-F1), intersection-based F1 score (IB-F1) are used to evaluate the model
performance. Median filters with fixed time length are used for post-
processing, and sound events have different thresholds from each other to
obtain hard predictions for calculating EB-F1. The metrics hyper-parameters
are identical to [12]. The model is trained using the Adam optimizer with a
maximum learning rate of 0.001, and ramp-up is used for the first 80 epochs.

Table 1: _SED performance comparison between models using different dynamic convolution on the validation set._ Model |  PSDS1 $\uparrow$ |  PSDS2 $\uparrow$ |  CB-F1 $\uparrow$ |  IB-F1 $\uparrow$
---|---|---|---|---
Baseline [19] | 0.370 | 0.579 | 0.469 | 0.714
DDFConv [18] | 0.387 | 0.624 | 0.467 | 0.720
FTDConv | 0.395 | 0.651 | 0.495 | 0.740
FFDConv | 0.436 | 0.685 | 0.526 | 0.751

###  3.2 Full-frequency dynamic convolution on SED

We compared the performances of baseline with full dynamic convolution
methods, including decoupled dynamic convolution (DDFConv) [18], full-temporal
dynamic convolution (FTDConv), and full-frequency dynamic convolution
(FFDConv). For full dynamic convolution methods, dynamic convolution layers
replaced all convolution layers except the first layer from the baseline model
[19].

The results are shown in Table 1. Three types of full dynamic convolution can
all outperform the baseline, which proves full dynamic convolution qualifies
in SED. In addition, it can be seen that the effects of three types of
convolution are in increasing order. First, FTDConv and FFDConv employ
content-adaptive temporal or frequency kernels, which can be viewed as giving
prior knowledge to SED compared with DDFConv. Second, FFDConv outperforms
FTDConv, which can prove that time-frequency patterns of sound events are
highly frequency-dependent, and this dependency is time-invariant. Moreover,
FFDConv models acoustic features with different kernels along the frequency
axis, which can be thought to be frequency components modeled in different
representation spaces. As if components of the feature are split into
different frequency spaces and then reassembled. This is consistent with the
characteristics of sound events.

Fig. 4: Feature comparison of FFDConv and CRNN. Features activation of the 5th
Conv block are shown in the 4th row. The trends of frequency band features
over time are shown in the 5th row. Note that y-axis labels of strong
prediction are abbreviations of the sound event categories. For example, Abr
stands for Alarm bell ringing.

###  3.3 Fine-grained modeling study

To explore FFDConv’s ability to understand acoustic spectral information at a
fine-grained level. We visualized feature of the middle layer. More
visualizations can be found in the supplementary material.

The visualization results are shown in Fig. 4. Comparing the features of
FFDConv and CRNN, we can see that most of the time-frequency patterns modeled
by CRNN are temporally isolated and disjoint. In contrast, FFDConv’s patterns
and their neighbors are in a whole, thereby forming a distinct time-frequency
representation. Moreover, this phenomenon can also be found in trends of
frequency band features over time. The waveforms of FFDConv are smoother than
CRNN. Specifically, the duration of peak and trough is longer in FFDConv’s
waveform, which results from the feature being mostly coherent over time.
There are more pulses in the resting state of CRNN’s waveforms, which are in a
disorganized state. Besides, the distributions of frequency band features are
consistent with alarm_bell_ring’s spectrogram in FFDConv’s waveforms. The
values of the low-frequency band features are smaller than those of the middle
and high-frequency bands when the alarm bell rings. However, the differences
between frequency bands in CRNN are ambiguous. As for the model’s prediction,
the CRNN’s isolated features directly lead to the incoherent output compared
with ground truth, which proves that the feature’s coherence over time is
essential. It’s interesting that the low-frequency white noise of the sound
clip is filtered by FFDConv, but CRNN tagged it as a speech. This has to do
with that dynamic convolution concentrates more on high-frequency texture
information, and white noise in the spectrogram lacks clear contour
information.

Actually, most SED models are trained in a frame-based supervised way, which
always leads to the feature and output being discrete over time. However,
FFDConv can alleviate this by frequency-dependent modeling, which models
different patterns for frequency bands, leading to a distinct representation
of a sound event. This modeling way is like an attention mechanism in which
the distribution of frequency band information of the spectrogram is
maintained. Besides, the convolution kernel for a frequency band is shared in
all frames, which produces temporally coherent representations. This is
consistent with both the continuity of the sound waveform and the vocal
continuity of sound events.

###  3.4 Ablation study

We compared the performance of different window sizes of the build kernel when
generating spatial filters. Note that the size of the spatial filter $K$ is
set to 3.

Table 2: _Comparison of different window size, W._ Model | $Atten$ | $W$ | PSDS1 | PSDS2
---|---|---|---|---
| ✗ | 3 | 0.421 | 0.650
| ✓ | 1 | 0.421 | 0.659
FFDConv | ✓ | 3 | 0.436 | 0.685
| ✓ | 5 | 0.423 | 0.656
| ✓ | 7 | 0.432 | 0.666

The results are shown in Table 2. With constraints of the attention module,
FFDConv can get better performance. This proves that before attention, spatial
filters of different frequency spaces may have a large distance from each
other. The performance of FFDConv is the best when window size is set to 3.
This is because the adjacent frequency components are considered compared to
size 1 when generating the spatial filter, and size 5 may suffer from
overfitting. In addition, it’s interesting that the performance recovers when
the window size is set to 7. This may have to do with the fact that dynamic
convolutions are relatively unstable.

##  4 conclusions

In this paper, we proposed full-frequency dynamic convolution, the first full-
dynamic method for SED. Full-frequency dynamic convolution is designed to
model time-frequency patterns in different frequency spaces. This design in
structure physically furnished 2D convolution with capability of frequency-
dependent modeling. Experiments on the DESED show that full-frequency dynamic
convolution is superior to not only baseline but also other full-dynamic
convolutions, which proves FFDConv qualifies in SED. In addition, by
visualizing features of sound events, we found that FFDConv can extract
temporally coherent features in specific frequency bands, which is consistent
with the vocal continuity of sound events. This proves that FFDConv has great
frequency-dependent perception ability. In the future, we aim to explore new
methods to model vocal continuity of sound events.

## References

  * [1] J. Rouat,  “Computational Auditory Scene Analysis: Principles, Algorithms, and Applications,”  IEEE TNN, vol. 19, no. 1, pp. 199–199, Jan. 2008.
  * [2] Zhengxia Zou, Keyan Chen, Zhenwei Shi, Yuhong Guo, and Jieping Ye,  “Object Detection in 20 Years: A Survey,”  Proc. IEEE, vol. 111, no. 3, pp. 257–276, March. 2023.
  * [3] Jiaxing Sun, Yujie Li, Huimin Lu, Tohru Kamiya, and Seiichi Serikawa,  “Deep Learning for Visual Segmentation: A Review,”  in COMPSAC, July. 2020, pp. 1256–1260.
  * [4] Qin Jin, Peter Schulam, Shourabh Rawat, Susanne Burger, Duo Ding, and Florian Metze,  “Event-based video retrieval using audio,”  in INTERSPEECH, Sept. 2012, pp. 2085–2088.
  * [5] Christian Debes, Andreas Merentitis, Sergey Sukhanov, Maria Niessen, Nikolaos Frangiadakis, and Alexander Bauer,  “Monitoring Activities of Daily Living in Smart Homes: Understanding human behavior,”  IEEE SPM, vol. 33, no. 2, pp. 81–94, March. 2016.
  * [6] Juan Pablo Bello, Charlie Mydlarz, and Justin Salamon,  Sound analysis in smart cities, pp. 373–397,  Springer International Publishing, Sept. 2017.
  * [7] Jie Hu, Li Shen, and Gang Sun,  “Squeeze-and-Excitation Networks,”  in CVPR, Jun. 2018, pp. 7132–7141.
  * [8] Xiang Li, Wenhai Wang, Xiaolin Hu, and Jian Yang,  “Selective Kernel Networks,”  in CVPR, Jun. 2019, pp. 510–519.
  * [9] Sanghyun Woo, Jongchan Park, Joon-Young Lee, and In So Kweon,  “CBAM: Convolutional Block Attention Module,”  in ECCV, Sept. 2018, pp. 3–19.
  * [10] Tong Na and Qinyi Zhang,  “Convolutional network with conformer for semi-supervised sound event detection,”  Tech. Rep., DCASE2021 Challenge, 2021.
  * [11] Xu Jia, Bert De Brabandere, Tinne Tuytelaars, and Luc V Gool,  “Dynamic filter networks,”  in NEURIPS, Dec. 2016.
  * [12] Hyeonuk Nam, Seong-Hu Kim, Byeong-Yun Ko, and Yong-Hwa Park,  “Frequency Dynamic Convolution: Frequency-Adaptive Pattern Recognition for Sound Event Detection,”  in INTERSPEECH, Sept. 2022, pp. 2763–2767.
  * [13] Chao Li, Aojun Zhou, and Anbang Yao,  “Omni-Dimensional Dynamic Convolution,”  in ICLR, Apr. 2022.
  * [14] Shengchang Xiao, Xueshuai Zhang, and Pengyuan Zhang,  “Multi-Dimensional Frequency Dynamic Convolution with Confident Mean Teacher for Sound Event Detection,”  in ICASSP, Jun. 2023, pp. 1–5.
  * [15] Julio Zamora Esquivel, Adan Cruz Vargas, Paulo Lopez Meyer, and Omesh Tickoo,  “Adaptive Convolutional Kernels,”  in ICCV Workshops, Oct. 2019, pp. 0–0.
  * [16] Zhi Tian, Chunhua Shen, and Hao Chen,  “Conditional Convolutions for Instance Segmentation,”  in ECCV, Aug. 2020, pp. 282–298.
  * [17] Xinlong Wang, Rufeng Zhang, Tao Kong, Lei Li, and Chunhua Shen,  “SOLOv2: Dynamic and Fast Instance Segmentation,”  in NEURIPS, Dec. 2020, pp. 17721–17732.
  * [18] Jingkai Zhou, Varun Jampani, Zhixiong Pi, Qiong Liu, and Ming-Hsuan Yang,  “Decoupled Dynamic Filter Networks,”  in CVPR, Jun. 2021, pp. 6647–6656.
  * [19] Emre Çakır, Giambattista Parascandolo, Toni Heittola, Heikki Huttunen, and Tuomas Virtanen,  “Convolutional Recurrent Neural Networks for Polyphonic Sound Event Detection,”  IEEE/ACM TASLP, vol. 25, no. 6, pp. 1291–1303, Jun. 2017.
  * [20] Antti Tarvainen and Harri Valpola,  “Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results,”  in NEURIPS, Dec. 2017.
  * [21] Hongyi Zhang, Moustapha Cisse, Yann N. Dauphin, and David Lopez-Paz,  “mixup: Beyond Empirical Risk Minimization,”  ICLR, Oct. 2018.
  * [22] Daniel S. Park, William Chan, Yu Zhang, Chung-Cheng Chiu, Barret Zoph, Ekin D. Cubuk, and Quoc V. Le,  “SpecAugment: A Simple Data Augmentation Method for Automatic Speech Recognition,”  in INTERSPEECH, Sept. 2019, pp. 2613–2617.
  * [23] Hyeonuk Nam, Seong-Hu Kim, and Yong-Hwa Park,  “Filteraugment: An Acoustic Environmental Data Augmentation Method,”  in ICASSP, May. 2022, pp. 4308–4312.


* Sandage & Tammann (1975) Sandage, A., & Tammann, G.A. 1975, ApJ, 197, 265
  * Santiago & Strauss (1992) Santiago, B.X. & Strauss, M.A. 1992, ApJ, 387, 9
  * Saslaw (2000) Saslaw, William C. 2000, The Distribution of the Galaxies: Gravitational Clustering in Cosmology, Cambridge University Press: Cambridge
  * Scargle (1998) Scargle, J. 1998, ApJ, 504, p.405.
  * Scargle (2002) Scargle, J. 2002 “Bayesian blocks in two or more dimensions: Image segmentation and cluster analysis,” in Bayesian Inference and Maximum Entropy Methods in Science and Engineering, American Institute of Physics Conference Proceedings, Volume 617, pp. 163-173.
  * Scargle et al. (2008) Scargle, J., Norris, J., & Jackson, B. 2008 “Studies in Astronomical Time Series Analysis. VI. Optimal Segmentation: Blocks, Triggers, and Histograms,” in preparation
  * Schaap & van de Weygaert (2000) Schaap, W.E. & van de Weygaert, R. 2000, A&A, 363, L29
  * Schaap (2007) Schaap, W.E. 2007, “The Delaunay Tessellation Field Estimator”, Ph.D. Thesis, Groningen University
  * Schlegel et al. (2009) Schlegel, D.J. et al. 2009, arXiv:0904.0468v3
  * Scott (1992) Scott, David W. 1992, Multivariate Density Estimation: Theory, Practice and Visualization, John Wiley & Sons, Inc.: New York
  * Shandarin (1983) Shandarin, S.F. 1983, Sov Astr. Letters, 9, 104
  * Shane & Wirtanen (1967) Shane, C.D. & Wirtanen, C.A. 1967, Publ. Lick. Obs. 22, Part 1
  * Shapley (1933) Shapley, H. 1933, PNAS, 19, 389
  * Shectman et al. (1996) Shectman, S.A. et al. 1996, ApJ, 470, 172
  * Sheth et al. (2003) Sheth, J.V., Sahni, V., Shandarin, S.F. & Sathyaprakash, B.S. 2003, MNRAS, 343, 22
  * Sheth & Tormen (2004) Sheth, J.V. & Tormen, G. 2004, MNRAS, 350, 1385
  * Silverman (1986) Silverman, B.W. 1986, Density Estimation for Statistics and Data Analysis, (Chapman & Hall; reprinted in 1998 by CRC Press: Boca Raton)
  * Slezak et al. (1990) Slezak, E., Bijaoui, A. & Mars, G. 1990, A&A, 227, 301
  * Slezak et al. (1993) Slezak, E., de Lapparent, V. & Bijaoui, A. 1993, ApJ, 409, 517
  * Snyder (1991) Snyder, Donald L. & Miller, Michael I. 1991, Random Point Processes in Time and Space, 2nd edition, Springer-Verlag: New York
  * Sousbie, Colombi & Pichon (2009) Sousbie, T., Columbi, S.& Pichon C. 2009, MNRAS, 393, 457
  * Sousbie (2010) Sousbie, T. 2010, submitted to MNRAS, arXiv:1009.4105
  * Sousbie, Pichon & Kawahara (2010) Sousbie, T., Pichon, C. & Kawahara, H. 2010, submitted to MNRAS, arXiv:1009.4104
  * Stein (1997) Stein, M.L. 1997 in Feigelson E.D., Babu G.J., eds., Statistical Challenges in Modern Astronomy II. Springer-Verlag, New York, p.166
  * Stoyan et al. (1985) Stoyan, D., Kendall, W.S. & Mecke, J. 1995, Stochastic Geometry and Its Applications, ed. John Wiley & Sons, Chichester
  * Springel et al. (2005) Springel, V., et al. 2005, Nature, 435, 629
  * Strauss et al. (2002) Strauss, M. A., et al. 2002 AJ, 124, 1810
  * Stril et al. (2010) Stril, A., Cahn, R., & Linder E.V. 2010, MNRAS, 404, 239
  * Szapudi & Szalay (1998) Szapudi, I & Szalay, A.S. 1998, ApJ, 494, 41
  * Totsuji & Kihara (1969) Totsuji, H. & Kihara, T. 1969, PASJ, 21, 221
  * Turner & Gott (1976) Turner, E.L. & Gott, J.R. 1976, ApJS, 32, 409
  * Turner et al. (1979) Turner, E.L., Aarseth, S.J., Gott, J.R., Blanchard, N.T. & Mathieu, R.D. 1979, ApJ, 228, 684
  * Ueda & Itoh (1997) Ueda, H. & Itoh, M. 1997, PASJ, 49, 131
  * van de Weygaert (1994) van de Weygaert R. 1994, A&A, 283, 361
  * van de Weygaert (2003) van de Weygaert R. 2003, “The Cosmic Foam: Stochastic Geometry and Spatial Clustering across the Universe,” Invited contribution in Proceedings it Statistical Challenges in Modern Astronomy III, eds. E.D. Feigelson & G.J. Babu, Springer-Verlag, pp. 175-196
  * van de Weygaert & Schaap (2009) van de Weygaert R. & Schaap, W. 2009, “The Cosmic Web: Geometric Analysis”, in Data Analysis in Cosmology, Lecture Notes in Physics, vol. 665, Eds V.J. Martinez, E. Saar, E. Mart nez-Gonz lez, and M.-J. Pons-Border a. Berlin: Springer, 2009., p.291-413
  * van de Weygaert & Aragón-Calvo (2009) van de Weygaert R. & Aragón-Calvo, M. 2003, “Geometry and Morphology of the Cosmic Web: Analyzing Spatial Patterns in the Universe,” Invited review ISVD09 (International Symposium on Voronoi Diagrams and Engineering), Copenhagen, Denmark. IEEE CPS, E3781, ed. F. Anton.
  * Villmann et al. (1997) Villmann, T., Der, R., Herrmann, M., & Martinetz, T., 1997, “Topology Preservation in Self-Organizing Feature Maps: Exact Definition and Measurement.” _IEEE Transactions on Neural Networks_ , 8, pp. 256-266.
  * Vogeley et al. (1994) Vogeley, M.S., et al. 1994, ApJ, 420, 525
  * Wright (1750) Wright, T. 1750, “An Original Theory or New Hypothesis of the Universe” New York: Elsevier 1971
  * York et al. (2000) York, D.G. et al. 2000, AJ, 120, 1579
  * Yoshioka & Ikeuchi (1989) Yoshioka, S. & Ikeuchi, S. 1989, ApJ, 341, 16
  * Zehavi et al. (2002) Zehavi, I. et al. 2002, ApJ, 571, 172
  * Zehavi et al. (2010) Zehavi, I. et al. 2010, arXiv:1005.2413
  * Zeldovich (1970) Zel‘dovich, Ya. B. 1970, A&A5, 84
  * Zeldovich et al. (1982) Zel‘dovich, Ya. B., Einasto, J., & Shandarin, S.F. 1982, Nature, 300, 407
  * Zhang et al. (2010) Zhang, Y., Springel, V. & Yang, X. 2010, arXiv:1006.3768
  * Zwicky (1957) Zwicky, F. 1957, PASP, 69, 518
  * Zwicky et al. (1961) Zwicky, F., Wield, P., Herzog, E., Karpowicz, M., & Kowal, C.T. 1961-68, Catalogue of Galaxies and Clusters of Galaxies, 6 volumes. Pasadena, California Institute of Technology


Tissue evolution: Mechanical interplay of adhesion, pressure, and
heterogeneity

Tobias Büscher1, Nirmalendu Ganai1,2, Gerhard Gompper1, Jens Elgeti1*,

1 Theoretical Soft Matter and Biophysics, Institute of Complex Systems and
Institute for Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich,
Germany

2 Department of Physics, Nabadwip Vidyasagar College, Nabadwip, Nadia 741302,
India

*<EMAIL_ADDRESS>
## Abstract

The evolution of various competing cell types in tissues, and the resulting
persistent tissue population, is studied numerically and analytically in a
particle-based model of active tissues. Mutations change the properties of
cells in various ways, including their mechanical properties. Each mutation
results in an advantage or disadvantage to grow in the competition between
different cell types. While changes in signaling processes and biochemistry
play an important role, we focus on changes in the mechanical properties by
studying the result of variation of growth force and adhesive cross-
interactions between cell types. For independent mutations of growth force and
adhesion strength, the tissue evolves towards cell types with high growth
force and low internal adhesion strength, as both increase the homeostatic
pressure. Motivated by biological evidence, we postulate a coupling between
both parameters, such that an increased growth force comes at the cost of a
higher internal adhesion strength or vice versa. This tradeoff controls the
evolution of the tissue, ranging from unidirectional evolution to very
heterogeneous and dynamic populations. The special case of two competing cell
types reveals three distinct parameter regimes: Two in which one cell type
outcompetes the other, and one in which both cell types coexist in a highly
mixed state. Interestingly, a single mutated cell alone suffices to reach the
mixed state, while a finite mutation rate affects the results only weakly.
Finally, the coupling between changes in growth force and adhesion strength
reveals a mechanical explanation for the evolution towards intra-tumor
heterogeneity, in which multiple species coexist even under a constant
evolutianary pressure.

## Introduction

Mutations change the cell fitness and thus its chance to survive and
proliferate [1]. Advantageous mutations are more likely to persist due to
natural selection, which drives the evolution of a tissue towards fitter cells
[2]. Cancer represents an example of evolution on a short time scale [3].
Furthermore, cancer is a multistep process, i.e. several mutations are needed
for a tumor in order to develop and become malignant [4]. Hence, tumorigenesis
might be expected to happen in a serial manner, i.e. a cell acquiring a
”beneficial” mutation and taking over the whole tissue. After some time, a
daughter cell acquires another mutation and again takes over. Interestingly,
however, tumors do not consist of a single cell type, but instead several
subpopulations coexist within the same tumor. This is called intra-tumor
heterogeneity [5].

Each mutation changes certain biochemical properties of a cell. This ranges
from misfunction in the error correction machinery during DNA replication and
disruptions in signaling pathways to epigenetic changes in the expression
level of certain proteins [1, 6, 7]. All these changes can also affect the
mechanical properties of the mutated cell, e.g. mutated cells which express
less adhesion proteins might be able to detach from the primary tumor more
easily [8], necessary to form metastases. On the other hand, mechanics feeds
back onto growth in several ways, e.g. increased apoptosis rate due to
mechanical stresses [9, 10] or dependence of the growth of tissue spheroids on
the properties of the surrounding medium [11, 12, 13].

It is the mechanical contribution to tissue development that we want to focus
on in this work. For mechanically regulated growth, homeostatic pressure plays
an important role [14]. In the homeostatic state, when apoptosis and division
balance each other, a tissue exerts a certain pressure onto its surrounding,
the homeostatic pressure $P_{\text{H}}$. The tissue is able to grow as long as
the external pressure $P$ is smaller than $P_{\text{H}}$. For the competition
between different tissues for space, it has been suggested that the tissue
with the higher homeostatic pressure grows at the expense of the weaker
tissue. Several theoretical studies employ this concept in order to describe
interface propagation between two competing tissues [15, 16, 17]. A metastasis
would need to reach a critical size, below which the additional Laplace
pressure due to surface tension would cause the metastasis to shrink and
disappear [14]. However, reduced adhesion between tissues, which increases
surface tension, leads to an enhanced growth rate at the interface between
them, stabilizing coexistence even for differing homeostatic pressures [18].

In this work, we study the influence of mutations that change the mechanical
properties of cells on the competition dynamics, especially the interplay
between changes in the adhesive properties and the strength with which a cell
pushes onto its surrounding. Particularly interesting is the case where loss
of adhesion comes at the cost of lower growth strength. This is motivated by
the observed down-regulation of E-cadherin, an adhesion protein in epithelia,
in many types of cancer [19]. Interestingly, E-cadherin is also involved in
signaling processes connected to cell growth [20]. We find that in this case
several cell types with different mechanical properties can coexist and that
the cell type with the highest homeostatic pressure does not necessarily
dominate the competition.

## Results

Several models have been developed previously in order to study tissue growth
[21], in combination with different simulation techniques, including vertex
[22, 23] and particle-based [24, 25] models as well as Cellular Potts models
[26, 27]. We employ the two particle growth (2PG) model of Refs. [28, 29, 18].
A cell is described by two particles which repel each other via a growth force

$\textbf{{F}}_{ij}^{\text{G}}=\frac{G}{(r_{ij}+r_{0})^{2}}\hat{\textbf{{r}}}_{ij}\text{,}$
(1)

with strength $G$, unit vector $\hat{\textbf{{r}}}_{ij}$, distance $r_{ij}$
between the two particles and a constant $r_{0}$. Different cells interact via
a soft repulsive force $\textbf{{F}}_{ij}^{\text{V}}$ on short distances,
maintaining an excluded volume, and a constant attractive force
$\textbf{{F}}_{ij}^{\text{A}}$ on intermediate distances, modeling cell-cell
adhesion, with

$\left.\begin{array}[]{@{}ll@{}}\textbf{{F}}_{ij}^{\text{V}}&=f_{0}\left(\frac{R_{\text{PP}}^{5}}{r_{ij}^{5}}-1\right)\hat{\textbf{{r}}}_{ij}\\\
\textbf{{F}}_{ij}^{\text{A}}&=-f_{1}\hat{\textbf{{r}}}_{ij}\par\end{array}\right\\}\text{for
}r_{ij}<R_{\text{PP}}\text{,}$ (2)

with exclusion coefficient $f_{0}$, adhesion strength coefficient $f_{1}$, and
cut-off length $R_{\text{PP}}$. A cell divides when the distance between its
two particles reaches a size threshold $r_{\text{ct}}$. A new particle is then
placed close (randomly within a short distance $r_{\text{d}}$) to each of the
two particles of the divided cell. Each of these pairs then constitutes a new
cell. Apoptosis is modeled by removing cells randomly at a constant rate
$k_{\text{a}}$.

We employ a dissipative particle dynamics-type thermostat, with an effective
temperature $T$, to account for energy dissipation and random fluctuations. We
choose the value of $T$ such that cells can escape local minima, but other
thermal effects are negligible. Note that all parameters can be set
individually for each cell type as well as between different cell types for
inter-cell interactions. We only vary the growth-force strength $G^{\alpha}$
and adhesion strength $f_{1}^{\alpha\beta}$ between cells of the same
($\alpha=\beta$) and different ($\alpha\neq\beta$) cell types, respectively,
where $\alpha$ and $\beta$ are cell-type numbers. We report simulation
parameters relative to a standard host cell type (see Materials and methods
for numerical values), denoted with a dagger, e.g. $G^{\dagger}=G/G^{0}$. Time
is measured in terms of the inverse apoptosis rate $k_{\text{a}}$, distance in
units of the pair potential cut-off length $R_{\text{PP}}$ and stresses in
units of $G^{0}/R_{\text{PP}}^{4}$. Quantities reported in these units are
denoted by an asterisk ∗. All simulations are performed in a cubic box with
edge length $L=12\cdot R_{\text{PP}}$ and periodic boundary conditions in all
directions, unless stated otherwise.

Fig 1: Evolution of a tissue with mutations altering growth-force strengh
$G^{\dagger}$ and adhesion strength $f_{1}^{\dagger}$ independently. Heatmaps
displaying cell-number fractions $\phi_{\alpha}$ after a) zero generations
(initial condition), b) 50 generations, c) 100 generations and d) 125
generations.

Tumor cells even within the same tumor are not all identical, but vary in
terms of all kind of attributes, e.g. expression levels of different proteins
[30] or their reaction to certain treatments [31]. Hence, there is not only a
competition between the tumor and the host, but also between cell-
subpopulations of the tumor. Different models exist to describe tumor
heterogeneity, e.g. cancer stem cells [32] or clonal evolution [33]. In the
latter case, a tumor originates from a single mutated cell, which can acquire
additional mutations over time, yielding additional subpopulations. We model
this behaviour by defining a fixed number $n$ of different ”genotypes”, each
having a different growth-force strength $G^{\alpha}$ and adhesion strength
$f_{1}^{\alpha\alpha}$. Mutations are implemented by offering each daughter
cell after a division event the chance to change its genotype with a certain
probability.

In tissues, several adhesion mechanisms exist, serving a variety of different
functions to maintain tissue integrity. Between epithelial cells, the strength
of cell-cell adhesion is to a large degree regulated by anchoring junctions,
e.g. adherens junctions, which connect the actin cytosceletons of neighbouring
cells. Adherens junctions are mediated by cadherins, which form homophilic
bonds between cells. Thus, the strength of adhesion between cells is limited
by the cell expressing less cadherin, or, in terms of our simulation model
$f_{1}^{\alpha\beta}=\min(f_{1}^{\alpha\alpha},f_{1}^{\beta\beta})$. A reduced
adhesion strength yields a higher homeostatic pressure [29], which is
otherwise dominated by the growth-force strength $G$. For free parameter
evolution, the tissue thus evolves to a strong-growing and low-adhesive
genotype (see Fig 1), as predicted by the homeostatic pressure approach [14].

Fig 2:  Time evolution of the cell-number fractions $\phi_{\alpha}$ of each
genotype for tradeoff paramter a) $\tau=0$, b) $\tau=1$, c) $\tau=2$ and d)
$\tau\rightarrow\infty$, $d\rightarrow 0$. Simulations start from a host
(standard) tissue at homeostasis, with $n=21$ genotypes, $p_{\text{m}}=0.01$
in all and $d=0.025$ in a)-c). White space corresponds to times where no cells
of the genotype exist. Color is coded on a logarithmic scale. Curves above
display homeostatic pressure $P_{\text{H}}^{\alpha*}$ (black solid), growth-
force strength $G^{\alpha\dagger}$(red dashed) and self adhesion strength
$f_{1}^{\alpha\alpha\dagger}$ (green dotted) of the corresponding genotype.

However, E-cadherin also plays a role in signaling processes connected to cell
growth, and thus a reduced expression might come at the cost of a lower
growth-force strength $G$, which in turn yields a lower homeostatic pressure.
We thus turn our attention to the case where an increase in growth-force
strength $G^{\alpha}$ comes at the cost of a higher self-adhesion strength
$f_{1}^{\alpha\alpha}$. We assume the relations as

$\displaystyle G^{\alpha}$ $\displaystyle=(1+D^{\alpha})G^{0}$ (3)
$\displaystyle f_{1}^{\alpha\alpha}$
$\displaystyle=(1+D^{\alpha}\cdot\tau)f_{1}^{0}\text{,}$ (4)

with genotype number $\alpha$ in the range $[-(n-1)/2,(n-1)/2]$, evolutionary
distance $D^{\alpha}=d\cdot\alpha$, distance $d$ between neighbouring
genotypes and tradeoff paramteter $\tau$ (with
$G^{\alpha},f_{1}^{\alpha\alpha}>0\ \forall\ \alpha$). After a division event,
each daughter cell might mutate into a new genotype with probability
$p_{\text{m}}$. If the cell mutates, its genotype number is changed to
$\alpha_{\text{mother}}~{}\pm 1$ randomly. This yields a mutation rate
$k_{\text{m}}=2p_{\text{m}}k_{\text{a}}$.

Figure 2 displays results of such simulations for four different cases: only
variation of growth-force strength ($\tau=0$), balanced tradeoff ($\tau=1$),
adhesion strength varied twice as much as growth-force strength ($\tau=2$) and
only variation of adhesion strength ($\tau\rightarrow\infty$). Without
tradeoff (Fig 2a)), the tissue evolves towards the strongest growing genotype
or, equivalently, the one with the highest homeostatic pressure. Similarly,
for $\tau\rightarrow\infty$ (Fig 2d)), the system evolves towards the lowest
adhesive genotype (again, the one with the highest $P_{\text{H}}$). We find
the most dynamic evolution for a balanced tradeoff (Fig LABEL:snapshot and
2b)). At first, the system evolves to stronger growing and more adhesive
genotypes. Over time a noticable fraction of cells evolves also towards weak-
growing, less adhesive genotypes. The cell-number fractions
$\phi_{\alpha}=N_{\alpha}/N$ (with individual and total number of cells,
$N_{\alpha}$ and $N$), show large fluctuations (see Fig LABEL:snapshotb) and
c)), with individual genotypes not being populated at all for certain time
periods. Besides this highly dynamic temporal evolution, after an initial time
period the system is dominated by genotypes with increased growth force and
adhesion strength at all times, with the one at the upper boundary having the
highest cell-number fraction for most of the time (see Fig LABEL:snapshota)).
This result comes at a surprise, as this is also the genotype with the lowest
homeostatic pressure, while the one at the lower boundary, which is basically
never populated, has the highest $P_{\text{H}}$. For a higher tradeoff (Fig
2c)), we still find a broad distribution of genotypes, with less adhesive
genotypes dominating over the stronger growing ones, i.e. the loss in growth-
force strength is overcompensated by a lower adhesion strength.

In order to gain insight into the underlying mechanism of this dynamic
evolution, we study the competition between two genotypes and no mutations
($p_{\text{m}}=0$). Simulations are started from a single mutated cell (with
increased/decreased growth force and adhesion strength) in a host tissue at
the homeostatic state (we label the mutant with M and the host (wild type)
with W). Even in this simplified case, we find one parameter regime in which
the mutant is not able to grow, one regime with stable coexistence in a highly
mixed state and another regime in which the mutant outcompetes the host.
Figure LABEL:number_fractions shows the averaged number fractions of the
mutant at the steady state. For reduced growth force and adhesion strength
(Fig LABEL:number_fractionsa)), the mutant can only grow against the host if
its adhesion strength is reduced below a critical $f_{1}^{\text{crit}}$. In
terms of Eq. (4), the value of $f_{1}^{\text{crit}}$ roughly corresponds to a
balanced tradeoff ($\tau\approx 1$). Already for
$f_{1}^{\text{MM}}>f_{1}^{\text{crit}}$, the homeostatic pressure of the
mutant exceeds the one of the host, i.e. a parameter regime exists in which
the mutant is not able to grow, despite of the higher $P_{\text{H}}$. The
reverse happens when growth force and adhesion strength are increased. The
mutant completely takes over the compartment, although its homeostatic
pressure is smaller than that of the host. Again, coexistence is only found
when the adhesion strength is increased above $f_{1}^{\text{crit}}$. In the
coexistence regime, the mutant number fraction scales as
$\phi^{\text{M}}\propto 1/(f_{1}^{\text{MM}}-f_{1}^{\text{WW}})$.

Altogether, the competition between two genotypes alone yields the same
qualitative results as the more complex multi-genotype case discussed before.
Still, the question remains how a genotype with lower homeostatic pressure can
outcompete a stronger genotype. The answer can only lie in the adhesion
strength $f_{1}^{\text{MW}}=\min(f_{1}^{\text{MM}},f_{1}^{\text{WW}})$ between
mutant and host cells. This choice of cross-adhesion strength breaks symmetry,
as the stronger adhering genotype has more free space at the interface, which
favors divisions [18].

To address this question, we develop a phenomenological model which
incorporates pressure-dependent growth as well as interfacial effects, in
order to obtain a qualitative explanation of the simulation results.

We start with the expansion of the bulk growth rate $k_{\text{b}}$ around the
homeostatic pressure,

$k_{\text{b}}=\kappa(P-P_{\text{H}})\text{,}$ (5)

with the pressure response coefficient $\kappa$. Due to the high degree of
mixing, the number fractions $\phi^{\text{M/W}}$ and hence the strengh of
interfacial effects vary locally. In a mean-field approximation, we take the
interfacial effects to be proportional to
$\phi^{\text{M}}(1-\phi^{\text{M}})$, with individual prefactors $\Delta
k_{\text{s}}^{\text{M/W}}$ for each genotype. The time evolution is then given
by

$\displaystyle\partial_{t}\phi^{\text{M}}=$
$\displaystyle\kappa(P_{\text{H}}^{\text{M}}-P)\phi^{\text{M}}+\Delta
k_{\text{s}}^{\text{M}}\phi^{\text{M}}(1-\phi^{\text{M}})$ (6)
$\displaystyle\partial_{t}(1-\phi^{\text{M}})=$
$\displaystyle\kappa(P_{\text{H}}^{\text{M}}+\Delta
P_{\text{H}}-P)(1-\phi)+\Delta
k_{\text{s}}^{\text{W}}\phi^{\text{M}}(1-\phi^{\text{M}})\text{,}$ (7)

with the difference in homeostatic pressure $\Delta
P_{\text{H}}=P_{\text{H}}^{\text{W}}-P_{\text{H}}^{\text{M}}$. Addition of
Eqs. (6) and (7) yields the pressure

$\displaystyle P=P_{\text{H}}^{\text{W}}-\Delta
P_{\text{H}}\phi^{\text{M}}+\frac{\Delta k_{\text{s}}^{\text{M}}+\Delta
k_{\text{s}}^{\text{W}}}{\kappa}\phi^{\text{M}}(1-\phi^{\text{M}})\text{.}$
(8)

Thus, the pressure is given by the homeostatic pressures of the two genotypes
weighted by their number fraction plus an interfacial term. A figure
displaying the pressure measured during the simulations shown in Fig
LABEL:number_fractions can be found in the S1 Appendix.. Insertion of Eq. (8)
into Eq. (6) yields a differential equation for the number fraction with three
fixed points ($\partial_{t}\phi^{\text{M}}=0$), $\phi_{1}^{\text{M}}=0$,
$\phi_{2}^{\text{M}}=1$, and

$\displaystyle\phi_{3}^{\text{M}}=\frac{-\kappa\Delta P_{\text{H}}+\Delta
k_{\text{s}}^{\text{M}}}{\Delta k_{\text{s}}^{\text{M}}+\Delta
k_{\text{s}}^{\text{W}}}\text{.}$ (9)

We discuss this result for the case of reduced growth force and adhesion
strength of the mutant. $\Delta k_{\text{s}}^{\text{M}}$ might be expected to
vanish, as $f_{1}^{\text{MM}}=f_{1}^{\text{MW}}$ and mutant cells thus would
not feel whether neighbouring cells are mutant or host cells. However, in
order to grow, a cell needs to impose a strain on its surrounding. Host cells
adhere more strongly to each other, thus it is harder for a mutant cell to
impose a strain when surrounded by host cells. Hence, $\Delta
k_{\text{s}}^{\text{M}}$ is actually negative and the homeostatic pressure of
the mutant needs to exceed the host pressure by $-\Delta
k_{\text{s}}^{\text{M}}/\kappa$ in order to be able to grow against the host.
At this point, $\phi_{3}^{\text{M}}$ becomes positive, as long as $\Delta
k_{\text{s}}^{\text{M}}+\Delta k_{\text{s}}^{\text{W}}>0$. Host cells can
impose a strain more easily when surrounded by mutant cells and, additionally,
have more free space than when surrounded by other host cells. Hence, $|\Delta
k_{\text{s}}^{\text{M}}|<\Delta k_{\text{s}}^{\text{W}}$ and the above
mentioned condition is fulfilled. Similarly, coexistence can be found for
increased growth force and adhesion strength when $\Delta P_{\text{H}}>-\Delta
k_{\text{s}}^{\text{W}}/\kappa$. The above mentioned scaling of the mutant
number fraction can be obtained by an expansion of $\Delta P_{\text{H}}$ and
$\Delta k_{\text{s}}^{\text{M/W}}$ to linear order in terms of
$\epsilon:=(f_{1}^{\text{MM}}-f_{1}^{\text{WW}})/f_{1}^{\text{WW}}$ in Eq.
(9),

$\displaystyle\phi_{3}^{\text{M}}=\frac{-\kappa\Delta
P_{\text{H}}^{0}}{(\Delta k_{\text{s}}^{\text{M1}}+\Delta
k_{\text{s}}^{\text{W1}})\epsilon}+\frac{-\kappa\Delta P_{\text{H}}^{1}+\Delta
k_{\text{s}}^{\text{M1}}}{\Delta k_{\text{s}}^{\text{M1}}+\Delta
k_{\text{s}}^{\text{W1}}}\text{.}$ (10)

The zeroth order terms of $\Delta k_{\text{s}}^{\text{M/W}}$ vanish as there
are no interfacial effects when the adhesion strength between host and mutant
cells is equal to their self-adhesion strength, while $\Delta
P_{\text{H}}^{0}$ can be non-zero due to a changed growth-force strength.
Indeed, Eq. (10) reproduces the simulation data reasonably well (see Fig
LABEL:number_fractions). A discussion of the numerical values of the fitted
parameters and additional results can be found in S1 Appendix..

Figure LABEL:slope displays similar results as shown in Fig
LABEL:number_fractions, but now as a function of the tradeoff $\tau$ in Eq.
(4). For $\tau<1$ the genotype with higher growth-force strength outcompetes
the weaker genotype, for $1<\tau<2$ a transition towards the less adhesive
genotype occurs, while for even higher values of the tradeoff $\tau>2$ the
less adhesive genotype outcompetes the second genotype. This transition from
strongly growing, adhesive to weakly growing, less adhesive genotypes is found
in the same range of $\tau$ as in the competition between many genotypes.
Hence, the simplified case of two competing genotypes captures the essential
physics to explain the coexistence between many competing genotypes and,
additionally, provides a quantitative description.

Next, we turn our attention to the effect of a finite mutation rate on the
evolution of the system. Figure LABEL:mutationrate_clustera) shows the number
fraction of the mutant as a function of $k_{\text{m}}$ for different
combinations of evolutionary distance $D^{\alpha}$ and tradeoff $\tau$, in
comparision to the number fraction reached for a single mutation event. As
expexted, the number fraction converges towards $1/2$ with increasing
$k_{\text{m}}$ for all combinations. For moderate mutation rates, however, the
number fraction largely fluctuates around the same average as of a single
mutation event. The single mutation leads to a stable coexistence of the two
genotypes - additional mutations quickly relax back to this state.
Siginificant deviations occur only if in the steady state of the single
mutation event the number fraction of the weaker genotype is close to zero. In
that case, the weaker genotype consists only of one or very few small cohesive
clusters of cells, because cells of the weaker genotype need to detach from
the primary cluster in order to form new clusters, but are likely to die when
they do so, as they are only surrounded by cells of the stronger genotype.
Hence, the distribution of cells is highly non-homogenous. Compared to the
single mutation event, even a small mutation rate leads to the formation of
multiple small cluster all over the system, thus increasing the number
fraction of the weaker genotype (see Fig LABEL:mutationrate_clusterb) for
comparision in terms of number of clusters and Materials and methods for
further discussion). This result explains why at least two genotypes, in
addition to the dominating genotype, are populated as well in the cases shown
in Fig 2a) and d). When the number fractions of both genotypes are
sufficiently large (for $1\leq\tau\leq 2$), deviations from the average of a
single mutation are still small for the standard mutation probability.
Additionally, in the competitions between many genotypes, mutations change the
genotype to $\alpha\pm 1$ randomly and not in a preferred direction. Hence, we
conclude that the precise value of the mutation probability does not play an
important role in the regime where we find a heterogeneous distribution of
genotypes, as long as it is reasonably small ($k_{\text{m}}\ll k_{\text{a}}$).

Given that a single mutated cell can grow to tissue of macroscopic size in a
certain parameter regime for
$f_{1}^{\text{MW}}=\min(f_{1}^{\text{MM}},f_{1}^{\text{WW}})$, the question
arises how likely it is to actually reach this state. In order to study this
probability, we mutate again a single cell in a host tissue at its homeostatic
state. A mutation that reaches a certain threshold $N_{\text{t}}=20$ of cells
counts as a survival event (the chance to die after reaching this treshold
becomes extremely small), apoptosis of the last mutant cell as a death event.
Figure LABEL:survival shows the averages of many such simulations. For reduced
growth force and adhesion strength, the survival probability $p_{\text{s}}$ is
only non-zero below the critical adhesion strength $f_{1}^{\text{crit}}$. For
$f_{1}^{\text{MM}}<f_{1}^{\text{crit}}$, $p_{\text{s}}$ increases linearly
with further decreasing adhesion strength. On the other hand, when growth
force and adhesion strength are increased, the survival probability first
shows a plateau, whose value increases with increasing growth force strength,
from which it will probably drop to zero with further increase. Simulations in
this regime are difficult, because a mutated cell can easily grow to a few
cells, but will hardly reach the number threshold nor completely vanish again.
Due to the high self-adhesion strength on the one hand, it becomes hard to
detach from the other cells, but on the other hand easy to grow against the
host when only few or no other mutant cells are around. This explains the
larger error bars at the highest values of the adhesion strength, where the
sample size is small.

## Discussion

We have shown how intra-tumor heterogeneity, the existence of multiple
subpopulations within the same tumor, can arise due to mechanical interactions
alone. The simultaneous change of the adhesion and growth-force strength
stabilizes the coexistence of multiple subpopulations, in a highly dynamic
state. A higher growth-force strength alone, as well as a lower adhesion
strength, favor proliferation of a single subpopulation and the evolution of
the system to cell types with the highest growth-force strength, or lowest
adhesion strength, respectively. A tradeoff between the two, however, yields
coexistence between multiple subpopulations of different cell types.
Interestingly, the expression of the adhesion protein E-cadherin, which also
affects cell growth, has been found to be down-regulated in many real tumors
[19].

The simulations also reveal that the homeostatic pressure of a cell type is
not necessarily the only quantity that determines the result of a competition.
Interactions between different cell types, in our model determined by the
adhesion between them, can lead to a completely reverse outcome, i.e. a cell
type with lower homeostatic pressure can outcompete a stronger one completely.
A phenomenological model explains the results on a qualitative level. The
evolution of each cell type is governed by mechanically-regulated growth,
while mutation rates only play a minor role in the dynamics.

An interesting future aspect to be studied is the influence of open
boundaries. A tissue with a negative homeostatic pressure then naturally grows
to a spheroid of finite size, with an enhanced rate of division at the surface
[29]. For competing cell types, this would lead to an interplay between
surface and interfacial effects.

## Materials and methods

### Standard (host) tissue and simulation parameters

We define a set of reference simulation parameters, which we refer to as host
parameters. Table 1 shows the values in simulation units. In simulations we
keep the host W fixed and vary the parameters of the mutant M around the
values of the host.

Table 1: Simulation parameters and measured properties of the standard (host)
tissue

.  Parameter Symbol Value Time Step $\Delta t$ $10^{-3}$ Pair potential
interaction range $R_{\text{PP}}$ $1$ Cellular expansion pressure constant
$r_{0}$ 1 Cell division distance treshold $r_{\text{ct}}$ $0.8$ New cell
particle initial distance $r_{\text{d}}$ 0.00001 Growth-force strength $G$
$40$ Mass $m$ 1 Intracell dissipation coefficient $\gamma_{\text{c}}$ $100$
Intercell dissipation coefficient $\gamma_{\text{t}}$ $50$ Background
dissipation coefficient $\gamma_{\text{b}}$ $0.1$ Apoptosis rate
$k_{\text{a}}$ 0.01 Mutation propability $p_{\text{m}}$ 0.01 Noise intensity
$k_{\text{B}}T$ 0.1 Repulsive cell-cell potential coefficient $f_{0}$
$2.39566$ Attractive cell-cell potential coefficient $f_{1}$ $6.0$ Isothermal
compressibility $\beta_{\text{T}}$ 1 Relaxation time constant $t_{\text{P}}$
$1$ Homeostatic pressure $P_{\text{H}}^{*}$ $0.1321\pm 0.0005\text{\,}$
Pressure response coefficient $\kappa^{*}$ $2.676\pm 0.080\text{\,}$

### Cluster analysis

As explained in the results section, a constant rate of mutation leads to an
enhanced formation of clusters when the weaker genotype is barely able to grow
against the stronger genotype and consists of only one or few clusters for a
single mutation event. We define a cluster as all cells of the same genotype
that are in interaction range to at least one other member of the cluster
(DBSCAN clustering algorithm with number of minimal points equal to one).
Figure LABEL:mutationrate_clusterb) displays the number of clusters of the
weaker genotype in the competitions displayed in Fig
LABEL:mutationrate_clustera), in comparison to the result of a single mutation
event. Indeed, when the number fraction of the weaker genotype is small for a
the single mutation event ($\tau=1$), we find significant deviations even for
small mutation rates. In this case, the number of clusters first strongly
increases with mutation rate, with roughly a tenfold increase at the peak. For
even higher mutation probability, the number of clusters decreases again, due
to merging of clusters, finally leading to percolation.

## Supporting information

#### S1 Appendix.

Additional Results. Appendix containing additional results and the
corresponding figures.

## Acknowledgments

The authors gratefully acknowledge the computing time granted through JARA-HPC
on the supercomputer JURECA [35] at Forschungszentrum Jülich.

## References

  * 1. Weinberg RA.  The biology of cancer.  Garland Publishing, Inc.; 2007.
  * 2. Greaves M, Maley CC.  Clonal evolution in cancer.  Nature 2012;481(7381):306.
  * 3. Bozic I, Reiter JG, Allen B, Antal T, Chatterjee K, Shah P, et al.  Evolutionary dynamics of cancer in response to targeted combination therapy.  eLife. 2013;2:e00747.  doi:10.7554/eLife.00747.
  * 4. Vogelstein B, Kinzler KW.  The multistep nature of cancer.  Trends Genet. 1993;9(4):138 – 141.  doi:https://doi.org/10.1016/0168-9525(93)90209-Z.
  * 5. Heppner GH.  Tumor Heterogeneity.  Cancer Res. 1984;44(6):2259–2265.
  * 6. Preston BD, Albertson TM, Herr AJ.  DNA replication fidelity and cancer.  Semin Cancer Biol. 2010;20(5):281 – 293.  doi:https://doi.org/10.1016/j.semcancer.2010.10.009.
  * 7. Schnekenburger M, Diederich M.  Epigenetics Offer New Horizons for Colorectal Cancer Prevention.  Curr Colorectal Cancer Rep. 2012;8(1):66–81.  doi:10.1007/s11888-011-0116-z.
  * 8. Petrova YI, Schecterson L, Gumbiner BM.  Roles for E-cadherin cell surface regulation in cancer.  Mol Biol Cell. 2016;27(21):3233–3244.  doi:10.1091/mbc.E16-01-0058.
  * 9. Wernig F, Xu Q.  Mechanical stress-induced apoptosis in the cardiovascular system.  Prog Biophys Mol Bio. 2002;78(2):105 – 137.  doi:https://doi.org/10.1016/S0079-6107(02)00008-1.
  * 10. Cheng G, Tse J, Jain RK, Munn LL.  Micro-environmental mechanical stress controls tumor spheroid size and morphology by suppressing proliferation and inducing apoptosis in cancer cells.  PLoS One. 2009;4(2):e4632.  doi:10.1371/journal.pone.0004632.
  * 11. Montel F, Delarue M, Elgeti J, Malaquin L, Basan M, Risler T, et al.  Stress Clamp Experiments on Multicellular Tumor Spheroids.  Phys Rev Lett. 2011;107(18):188102.  doi:10.1103/PhysRevLett.107.188102.
  * 12. Alessandri K, Sarangi BR, Gurchenkov VV, Sinha B, Kiesling TR, Fetler L, et al.  Cellular capsules as a tool for multicellular spheroid production and for investigating the mechanics of tumor progression in vitro.  Proc Natl Acad Sci USA. 2013;110(37):14843–14848.  doi:10.1073/pnas.1309482110.
  * 13. Helmlinger G, Netti PA, Lichtenbeld HC, Melder RJ, Jain RK.  Solid stress inhibits the growth of multicellular tumor spheroids.  Nat Biotechnol. 1997;15(8):778–783.  doi:10.1038/nbt0897-778.
  * 14. Basan M, Risler T, Joanny JF, Sastre-Garau X, Prost J.  Homeostatic competition drives tumor growth and metastasis nucleation.  HFSP J. 2009;3(4):265–272.  doi:10.2976/1.3086732.
  * 15. Williamson JJ, Salbreux G.  Stability and Roughness of Interfaces in Mechanically Regulated Tissues.  Phys Rev Lett. 2018;121:238102.  doi:10.1103/PhysRevLett.121.238102.
  * 16. Ranft J, Aliee M, Prost J, Jülicher F, Joanny JF.  Mechanically driven interface propagation in biological tissues.  New J Phys. 2014;16(3):035002.
  * 17. Podewitz N, Jülicher F, Gompper G, Elgeti J.  Interface dynamics of competing tissues.  New J Phys. 2016;18(8):083020.
  * 18. Ganai N, Büscher T, Gompper G, Elgeti J.  Mechanics of tissue competition: interfaces stabilize coexistence.  New J Phys. 2019;21(6):063017.  doi:10.1088/1367-2630/ab2475.
  * 19. Beavon IRG.  The E-cadherin–catenin complex in tumour metastasis: structure, function and regulation.  Eur J Cancer. 2000;36(13):1607 – 1620.  doi:https://doi.org/10.1016/S0959-8049(00)00158-1.
  * 20. Pece S, Gutkind JS.  Signaling from E-cadherins to the MAPK Pathway by the Recruitment and Activation of Epidermal Growth Factor Receptors upon Cell-Cell Contact Formation.  J Biol Chem. 2000;275(52):41227–41233.  doi:10.1074/jbc.M006578200.
  * 21. Van Liedekerke P, Palm MM, Jagiella N, Drasdo D.  Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results.  Comp Part Mech. 2015;2(4):401–444.  doi:10.1007/s40571-015-0082-3.
  * 22. Farhadifar R, Röper JC, Aigouy B, Eaton S, Jülicher F.  The Influence of Cell Mechanics, Cell-Cell Interactions, and Proliferation on Epithelial Packing.  Curr Biol. 2007;17(24):2095 – 2104.  doi:https://doi.org/10.1016/j.cub.2007.11.049.
  * 23. Alt S, Ganguly P, Salbreux G.  Vertex models: From cell mechanics to tissue morphogenesis.  Phil Trans R Soc B. 2017;372:20150520.  doi:10.1098/rstb.2015.0520.
  * 24. Drasdo D, Höhme S.  A single-cell-based model of tumor growth in vitro: monolayers and spheroids.  Phys Biol. 2005;2(3):133–147.  doi:10.1088/1478-3975/2/3/001.
  * 25. Schaller G, Meyer-Hermann M.  Multicellular tumor spheroid in an off-lattice Voronoi-Delaunay cell model.  Phys Rev E. 2005;71(5 Pt 1):051910.
  * 26. Graner F, Glazier JA.  Simulation of biological cell sorting using a two-dimensional extended Potts model.  Phys Rev Lett. 1992;69:2013–2016.  doi:10.1103/PhysRevLett.69.2013.
  * 27. Szabó A, Merks RM.  Cellular Potts Modeling of Tumor Growth, Tumor Invasion, and Tumor Evolution.  Front Oncol. 2013;3:87.  doi:10.3389/fonc.2013.00087.
  * 28. Basan M, Prost J, Joanny JF, Elgeti J.  Dissipative particle dynamics simulations for biological tissues: rheology and competition.  Phys Biol. 2011;8(2):026014.  doi:10.1088/1478-3975/8/2/026014.
  * 29. Podewitz N, Delarue M, Elgeti J.  Tissue homeostasis: A tensile state.  EPL. 2015;109(5):58005.
  * 30. Marusyk A, Almendro V, Polyak K.  Intra-tumour heterogeneity: a looking glass for cancer?  Nat Rev Cancer. 2012;12:323–334.
  * 31. Marusyk A, Polyak K.  Tumor heterogeneity: Causes and consequences.  Biochim Biophys Acta. 2010;1805(1):105 – 117.  doi:https://doi.org/10.1016/j.bbcan.2009.11.002.
  * 32. Shackleton M, Quintana E, Fearon ER, Morrison SJ.  Heterogeneity in Cancer: Cancer Stem Cells versus Clonal Evolution.  Cell. 2009;138(5):822 – 829.  doi:https://doi.org/10.1016/j.cell.2009.08.017.
  * 33. Nowell P.  The clonal evolution of tumor cell populations.  Science. 1976;194(4260):23–28.  doi:10.1126/science.959840.
  * 34. Allen MP, Tildesley DJ.  Computer simulation of liquids.  Oxford: Clarendon Press; 1989.
  * 35. Jülich Supercomputing Centre.  JURECA: Modular supercomputer at Jülich Supercomputing Centre.  Journal of large-scale research facilities. 2018;4(A132).  doi:10.17815/jlsrf-4-121-1.