Self-Organization of Self-Clearing Beating Patterns in an Array of Locally Interacting Ciliated Cells Formulated as an Adaptive Boolean Network

TL;DR

This paper models the self-organization of ciliary beating patterns on multiciliated epithelia using adaptive Boolean networks, revealing how local interactions and mucus dynamics lead to realistic, robust patterns and patchy wave fields.

Contribution

It introduces a novel abstract model of cilia as adaptive Boolean networks driven by mucus interactions, explaining patchy wave patterns and modular organization in epithelia.

Findings

Deterministic update schemes produce more realistic dynamics.

Mucus distribution influences network topology and pattern formation.

Modularity explains patchy mucus wave patterns.

Abstract

The observed spatio-temporal ciliary beat patterns on multiciliated epithelia are suspected to be the result of self-organizing processes on various levels. Here, we present an abstract epithelium model at the pluricellular level, which intends to make the self-organization of ciliary beating patterns as well as of the associated fluid transport across the airway epithelium plausible. Ciliated cells are modeled in terms of locally interacting oscillating two-state actuators. The local interactions among these boolean actuators are triggered by seeded mucus lumps. In the course of a simulation the actuators' state and the associated mucus velocity field self-organize in tandem. We suggest to consider the dynamics on multiciliated epithelia in the context of adaptive (boolean) networks. Within the framework of adaptive boolean networks ciliated cells represent the nodes and as the mucus establishes the local interactions among nodes, its distribution determines the topology of the network. Furthermore, we present the results and insights from comprehensive parameter studies. The results show evidence that so called deterministic update schemes, which are meant to represent intercellular signaling, lead to more realistic and robust dynamics and may therefore be favored by nature. Finally, we suppose that unciliated cells introduce a modular network topology on ciliated epithelia causing the self-organization taking place simultaneously in each "ciliated module". This reasoning provides the first consistent explanation for the meaning of the observed patchy expression patterns of the mucus modulation wave fields. Modularity may therefore be seen as a modular construction plan of nature for ciliated epithelia whose number of cells range over several order of magnitudes.

Tables

Table 1. Table 1: The table shows how the functioning states spread across the different topologies and update schemes. The values represent the corresponding proportions of all functioning states in [%].

	DAU	RAU1	RAU2	SRAU1	SRAU2
USL	$84.2$	$0.0$	$0.0$	$0.0$	$0.0$
UHL	$1.2$	$0.0$	$0.0$	$0.0$	$0.0$
BHL	$0.5$	$0.0$	$0.0$	$0.0$	$0.0$
BHL+L	$3.7$	$1.6$	$0.0$	$0.9$	$7.8$

Table 2. Table 2: The table reports the paramater sets of stereotypical representatives, the typical appearance of attracting states, the transport speed and particularly, the direction of transport as well as of the mean wave propagation is visualized.

	Parameters	Appearance	$\left|{\overrightarrow{v}}_{g\infty }\right|$	${\overrightarrow{v}}_{g\infty }$	$\overrightarrow{k}$
USL - stereo1	(DAU,TO,64%,0.25,0%)	“frozen”	0.69
USL - stereo2	(DAU,TO,64%,0.25,0%)	“re-organizing”	0.5
UHL - stereo1	(DAU,VC,128%,1,0%)	“frozen”	0.38
UHL - stereo2	(DAU,TO,256%,1,0%)	“gliders”	0.13
BHL - stereo1	(DAU,TO,256%,1,0%)	poorly organized	0.2
BHL+L - stereo1	(RAU1,OP,256%,1,0%)	frozen crystal	0.5
BHL+L - stereo2	(DAU,OP,256%,1,0%)	frozen	0.5
BHL+L - stereo3	(DAU,TO,256%,1,0%)	modular, re-organizing	0.3

Table 3. Table 1: Parameter settings leading to highest transport rates, highest transport speeds and most structured states for each cell alignment.

	Update	Boundary	Mucus (%)	Energy $ϵ$	Unciliated (%)
	DAU	HC	256	0.25	0
USL	DAU	HC	4	0.25	0
	DAU	HC	128	0.75	0
UHL	DAU	VC	256	1.0	0
	DAU	OP	256	1.0	0
BHL	DAU	OP	32	1.0	0
	DAU	HC	256	1.0	0
	DAU	OP	256	1.0	0
BHL+L	SRAU2	HC	1	1.0	0
	RAU1	HC	64	1.0	0

Table 4. Table 2: Observable values of the best transporting, the fastest transporting and the best ordered state for each cell alignment.

	${\rho }_c$	$⟨\cos \theta ⟩$	$m\cdot \left|{\overrightarrow{v}}_{g\infty }\right|$	$\left|{\overrightarrow{v}}_{g\infty }\right|$	$\left|{\overrightarrow{v}}_{g0}\right|$
	24.0	0.95	$4\cdot {10}^3$	$0.69$	0.02
USL	21.7	0.85	${10}^2$	0.95	0.10
	24.2	0.95	$2\cdot {10}^3$	0.72	0.05
UHL	10.6	0.86	$3\cdot {10}^3$	0.46	0.10
	15.4	0.51	$2\cdot {10}^3$	0.27	0.02
BHL	5.9	0.14	$3\cdot {10}^2$	0.39	0.15
	16.3	0.57	$2\cdot {10}^3$	0.27	0.01
	11.7	0.93	$3\cdot {10}^3$	0.47	0.06
BHL+L	24.9	0.82	$10$	0.57	0.21
	34.0	0.98	${10}^3$	0.48	0.10

Equations

$${\psi }_i^{t+1}=f_i({\psi }_{j_1(i)}^t,{\psi }_{j_2(i)}^t...,{\psi }_{K_i(i)}^t).$$

$$\psi_{ij}^{t+1}=f_{ij}\big{(}\{\psi_{pq}^{t}\}_{pq\in n_{ij}},\{m_{pq}^{t}\}_{pq\in n_{ij}}\big{)}~{}.$$

$${\psi }_{ij}^{t+1}=\text{XOR}({\psi }_{ij}^t,c_{ij}^t).$$

$$c_{ij}^{t}=\text{NOR}(l_{0},l_{1})\cdot\Theta\big{[}u-\Delta w\big{]}=1~{}.$$

$$\Delta w=p\Delta V=k_{B}T\text{ln}\big{(}\frac{W_{n}}{W_{n-1}}\big{)}~{}.$$

$$\psi_{ij}^{t+1}=\text{XOR}\big{\{}\psi_{ij}^{t},\text{NOR}(l_{0},l_{1})\cdot\Theta\big{[}r-F_{U}(\Delta w)\big{]}\big{\}}.$$

$$F_U(u)=\frac{1}{\frac{ϵ}{1-ϵ}\exp (-u)+1},$$

$${\vec{v}}_{ij}(t)\doteq \frac{{\sum }_{p,q\in n_{ij}}({\vec{r}}_{ij}-{\vec{r}}_{pq})(m_{pq}^{t+1}-m_{pq}^t)}{{\sum }_{p,q\in n_{ij}}{m}_{pq}^t}$$

$${\vec{v}}_g(t^{\prime })\doteq \frac{{\sum }_{ij}{\vec{v}}_{ij}(t^{\prime }){\sum }_{p,q\in n_{ij}}{m}_{pq}^{t^{\prime }}}{{\sum }_{ij}{\sum }_{p,q\in n_{ij}}{m}_{pq}^{t^{\prime }}}.$$

$${\vec{v}}_g(t^{\prime })={\vec{v}}_{g0}+({\vec{v}}_g(t^{\prime })-{\vec{v}}_{g0})\doteq {\vec{v}}_{g0}+\Delta {\vec{v}}_g(t^{\prime }),$$

$$\cos {\theta }_{ij}\doteq \underset{t^{\prime }\to \infty }{\lim }\frac{⟨{\vec{v}}_{ij}(t^{\prime }){⟩}_{t^{\prime }}\cdot ⟨{\vec{v}}_g(t^{\prime }){⟩}_{t^{\prime }}}{|⟨{\vec{v}}_{ij}(t^{\prime }){⟩}_{t^{\prime }}|\cdot |⟨{\vec{v}}_g(t^{\prime }){⟩}_{t^{\prime }}|},$$

$${\rho }_c(t^{\prime })=\text{max}\left\{\frac{1}{2}\sum _{v=v_1...v_2}\text{C}(v,\gamma (v),t^{\prime })|\gamma (v)=a\cdot v+b;a,b,v\in \mathbb{R};0\le v_1,v_2\le I;0\le \gamma (v)\le J\right\}$$

$$P(k_j,k_i)=\sum _{f=0}^{f=f_N}|F(k_j,k_i,f){|}^2.$$

$$\text{C}(\Delta j,\Delta i,\Delta t)=\sum _{i,j,t^{\prime }}\Psi (i,j,t^{\prime })\cdot \Psi (i-\Delta i,j-\Delta j,t^{\prime }+\Delta t),$$

$${\vec{v}}_{wave}=\underset{\tau \to 0}{\lim }\frac{\vec{\sigma }(\tau )}{\tau },$$

Full text

I Introduction

Motivation for the present work is to contribute to the understanding of the fascinating and omnipresent phenomenon of mucociliary transport.

The epithelium of our airways constitutes a self-cleaning surface protecting our lungs from a variety of inhaled substances such as exhaust, dust, bacteria and other harmful substances of micro- and submicrometer size. These particles get entrapped by the mucus layer lining the inner surface of the tracheal and pulmonary airways, which is propelled by the coordinated oscillatory movement of millions of subjacent cilia. Cilia are hair-like protrusions of the cell membrane, an overview of their structure and function can be found in Linck (2009); Satir and Christensen (2007).

To gain insight into the mucociliary clearing mechanism, many experimental techniques on various length- and time-scales are conducted. Experimental studies typically focus either on a structural or functional component on a typical scale, like the remarkable insight into the detailed molecular structure of the axoneme (see e.g. Sui and Downing (2006); Burgess et al. (2003)), the orientation of the ciliary beating plane Satir and Christensen (2007), the distribution of the different epithelial cell types (e.g. et al. Oliveira (2003); Plopper et al. (1983)) or the observation of mucociliary phenomena on the pluricellular level Ryser et al. (2007). Even though much work has been done on various scales, our understanding of the basic mesoscopic function of the system still appears rather limited. Experimental methods face several challenges: until today, it is not possible to observe the details of the mucociliary dynamics in vivo. Further, it is difficult to observe the mucociliary clearing mechanism under controlled conditions. And finally, it is highly complex and laborious to simultaneously measure structural and functional parameters at different scales.

Mathematical models are perfectly suited to conduct parameter studies in order to determine the effects of structural and functional parameters on mucociliary phenomena and therefore, serve as an alternative approach for the investigation of the intriguing mucociliary phenomena. Considering the implementation of mucociliary interactions the existing models can roughly be divided into two classes.

One class of models prescribes the motion of cilia, including their coordination, and concentrates on the hydrodynamics and rheology of the system. In these models the action of cilia is modeled as distributed oscillating momentum source, such as oscillating envelope Ross and Corrsin (1974), traction layer Smith et al. (2007a), active porous medium or oscillating array of cilia represented by a distribution of hydrodynamic singularities along the cilia’s centerline Smith et al. (2007b). The aim of these studies is to elaborate the geometrical and rheological conditions under which the system achieves an efficient transport (e.g. Lee et al. (2011)).

More relevant to the context of the present study is a second class of models, aimed at the understanding of the emergence of the pattern of motion. In these models the details of the motion of cilia or their coordination are not prescribed explicitly. The coordinated behavior rather emerges in the course of self-organization, during which the system components interact locally, what drives the system from an initially uncoordinated state towards a globally coordinated state exhibiting cooperative behavior.

Self-organization may play a role on various time- and length-scales, at various levels for the generation of mucociliary transport. As discussed in Marshall (2010) the interactions between cilia-generated fluid flow and planar cell polarity (PCP) signaling may lead to self-organization during the morphogenesis, which establishes the alignment of the axonemes on the individual cells. Much of work has been done on the molecular level, concerning the molecular motors (dynein motor proteins), their organization in the axoneme as well as the hydrodynamics of the resulting model cilium (eg. Riedel-Kruse et al. (2007); Hilfinger and Jülicher (2008); Hilfinger et al. (2009)). The complex motion pattern of a single cilium is thought to result from the self-organization of many dynein motor proteins generating stresses on the elastic microtubules in the axoneme. The next higher level, the cellular level, concerns the formation of metachronal waves by the self-organized synchronization of cilia covering a single cell Elgeti and Gompper (2013). For a long time it remained an open question whether the synchronization is achieved by cellular signaling (membrane potentials and calcium waves), or if it emerges spontaneously, due to interactions between the individual cilia. Today, the computational models indicate that hydrodynamic coupling is sufficient to induce synchronization Elgeti and Gompper (2013); Mitran (2007); Gueron et al. (1997). It is conceivable, however, that calcium signaling is needed for fine tuning the synchronization Salathe (2007).

Here we propose a first step to the next higher level, the pluricellular level, on which the self-organization of ciliary activity among ciliated cells is thought to generate the global wave field and fluid transport on the airway epithelium. In order to make the self-organized beat patterns as well as self-organized fluid transport across the airway epithelium plausible, we present an oversimplified model, which is intended to represent a virtual self-cleaning epithelium. Ciliated cells are modeled as actuators alternating between two possible states representing ciliary oscillations. Interactions between the two-state actuators are mediated by discrete mucus lumps and enmeshed dust particles. Whenever it is possible, the mucus droplets get displaced by the action of actuators and they may block their motion in certain configurations, which is prescribed by local interaction rules. This highly simplified model based on locally interacting motors self-organizes towards a virtual self-cleaning epithelium: the initially randomly distributed phases erratically displace the mucus lumps at the beginning of the simulation. As time passes the motors self-organize, which is expressed by emergent global spatio-temporal structures, resembling the metachronal wavelets, which have been observed on the ciliated tracheal epithelium Ryser et al. (2007). These metachronal wavelets efficiently transport the mucus lumps into a well defined direction.

This paper is organized as follows. In Sec.II we introduce our abstract epithelium model based on symmetrically interacting two-state actuators. The local mucociliary interactions are first formulated in an intuitive way, before they are reformulated as logical update functions of the local state and mucus distribution. The goal of Sec.II is to formulate our epithelium model in terms of an adaptive boolean network. Within the framework of adaptive boolean networks nodes are represented by actuators and the topology of the network (in effect the links between the nodes) is determined by the mucus distribution. Furthermore, we introduce the involved model parameters encompassing in particular different morphologies, boundary conditions and update schemes. Thanks to the simplicity of our epithelium model, it is possible to investigate the impact of the model parameters on the system’s behaviour by conducting comprehensive parameter studies. The corresponding simulation data are presented in Sec.III. The simulation data have been screened for any suspicious coherences between the model parameters and the network dynamics as well as between the parameters and their corresponding attracting states. Furthermore, we aimed at finding the mechanisms driving the model efficiently towards properly self-cleaning states. Finally, in Sec.IV we discuss our findings. In particular, we hypothesize about the meaning of the discovered dynamical aspects (which seem to be universal for our specific model of locally interacting cells) for the real biological system.

II Model Description

II.1 Symmetrically Interacting Two-State Actuators

Our model is based on symmetrically interacting two-state actuators. Each actuator represents a ciliated epithelial cell. By arranging many actuators in a parquet-like manner, the model represents a ciliated epithelium. Cilia corresponding to the same cell are assumed to move synchronously back and forth. This back and forth motion of cilia bundles is incorporated by the alternation between the actuators’ two possible states, which is illustrated in Fig.1. The actual state of an actuator can thus be expressed by the boolean state variable $\psi\in\{0,1\}$.

II.2 Morphology of Virtual Epithelia

The morphology of the airway epithelium is thought to be the result of self-organizing processes during the morphogenesis Marshall (2010). An important characteristic of the morphology is the distribution of the orientation of the ciliary beating plane, which can be determined by the orientation of the microtubules in the axoneme Satir and Christensen (2007). The impact of the (dis-)orientation of the ciliary beating plane i.e. of the axonemal orientation on mucociliary transport has been discussed e.g. in Rutland and De Iongh (1993); De Iongh and Rutland (1989). It has been concluded that a disorganized ciliary orientation may be a primary cause for mucociliary dysfunction and vice versa. Studies attempting to quantify the axonemal orientation of respiratory cilia usually determine the ciliary beating plane on a few cells derived from nasal brushing (e.g. Rautiainen (1988)). To the best of our knowledge, possible differences between the intra- and intercellular orientation of the ciliary beating plane is lacking so far in the literature. Therefore, we first of all assume that cilia belonging to a cell have a coinciding beating plane. As it is conceivable that the ciliary orientation on a cell might be more strict than between cells, the actuators can be oriented vertically or horizontally, as illustrated in Fig.1. Note that the two different possible orientations are not meant being perpendicular to each other, but to represent two distinct orientations of the ciliary beating plane differing by an arbitrary angle.

Here, we consider three different conceivable parquet-like cell alignments shown in Fig.2. In the following we shall refer to these three cell alignments as (from the left to the right according to Fig.2) unidirected square lattice (USL), unidirected hexagonal lattice (UHL) and bidirected hexagonal lattice (BHL).

Another important characteristic of the morphology represents the population densities of ciliated and unciliated cells, whose role for the mucociliary dynamics has not been considered so far. According to electron microscopic studies Plopper et al. (1983); et al. Oliveira (2003) the area covered by ciliated cells roughly varies between one and two thirds of the total surface of the tracheal lining. Consequently, unciliated cells should be considered in pluricellular models. Therefore, we include the fraction of unciliated epithelial cells by the parameter $f\in[0,1]$, representing the proportion of randomly distributed empty sites in an array of actuators.

Finally, as indicated by the dashed grid lines in Fig.2 the virtual epithelium is represented by a two-dimensional array. For convenience, a site located at the $i$-th row and $j$-th column is denoted as $\psi_{ij}$. $\psi_{ij}$ either represents an actuator, then $\psi_{ij}\in\{0,1\}$, or indicates an empty site representing an unciliated cell, then $\psi_{ij}=NAN$. Accordingly, the state of an array having $I$ rows and $J$ columns can be denoted as: $\Psi=\{\psi_{ij}\}$, where $\psi_{ij}\in\{0,1,NAN\}$.

II.3 Boundary Conditions

Experiments aiming at the characterization of spatio-temporal features of mucociliary phenomena on the tracheal ciliary epithelium are based on the approach of excising a rectangular piece of the cylindric trachea (e.g. Ryser et al. (2007); Lee et al. (2005); Yi et al. (2002); Wong et al. (1993)), which changes the boundary conditions from cylindrical to open. The excision might influence the collective dynamical behavior of the system. Therefore, we consider four different boundary conditions, which are illustrated in Fig.3 and shall be referenced in the following as: open boundaries (OP), vertical cylindric boundaries (VC), horizontal cylindric boundaries (HC) and toric boundaries (TO).

II.4 Mucus

Mucus is discontinuously incorporated by randomly seeding mucus droplets of equal size. On the one hand, the discrete mucus model allows to simplify mucociliary interactions. On the other hand, it accounts for the fact that the mucus blanket is anyway made up of excreted mucus “flakes”, “plaques” or ”droplets“ Van As and Webster (1974), which may coalesce into a continuous layer, if their density is sufficiently large Geiser et al. (1997).

According to Fig.1 each actuator provides an empty field (0 or 1), which can be occupied by mucus droplets. On the other hand, empty sites (unciliated cells) provide two fields (0 and 1), which can be occupied by mucus droplets. Thus, the distribution of mucus droplets on the virtual epithelium is given by $M=\{m_{ij}\}$; if $\psi_{ij}=NAN$, then $m_{ij}=m_{0ij}+m_{1ij}$ .

II.5 System Update

In the course of a simulation the actuators are actuated sequentially. This means that only the state of the actuated actuator and its adjacent mucus configuration is updated while the states of adjacent actuators do not change. As soon as an actuator and its local mucus configuration has been updated, a subsequent actuator is updated for which the changes of the prior step are taken into account. We shall label the sequence of update steps by the “time” superscript $t$. Furthermore, $t^{\prime}\doteq t/N$ labels the update of the whole network, consisting of $N$ actuators, in the following. In the context of discrete dynamical systems such sequential update schemes are called asynchronous. An excellent overview of the impact of different update schemes on the network dynamics of random boolean networks is provided in Gershenson (2004a, 2002). In order to distinguish the different update schemes, we shall use a slightly modified form of the terms proposed by C. Gershenson Gershenson (2002) (since we are not dealing with random boolean networks). We shall use the following update schemes (consider the corresponding illustrations in Fig.4):

• •

**DAU - Deterministic Asynchronous Update

**This update mode corresponds to a pre-defined sequence, in which the nodes/cells in the lattice are addressed. Keep in mind, however, that determinism of addressing doesn’t mean determinism of actuators motion. Whether an actuator moves or not depends on the mucociliary interactions, which are on a high degree stochastic, as we shall see in the proceeding sections.

• •

**RAU - Random Asynchronous Update

**At each time step a single node is chosen randomly and updated. The random choice of nodes has been applied for sampling with and without replacement, as in Cornforth et al. (2005). Random sequences of nodes generated without and with replacement are denoted as RAU1 and RAU2, respectively.

• •

**SRAU - Semi-Random Asynchronous Update

**The update scheme SRAU uses a wavelike activation: at each time step a strip of actuators gets addressed. Inside of each strip the actuators get addressed according to the RAU1 scheme. As the choice of the strip is predefined but inside the strips the sequence of actuators is chosen randomly, this addressing scheme has a semi-random character. This wave-like activation has been implemented for a wave traveling from the left to the right, as well as for a wave traveling from the top to the bottom, which we denote as SRAU1 and SRAU2, respectively.

• •

**Mimicking Purcell’s Two-Hinged Low-Reynolds-Number-Swimmer

**As non-reciprocal motions play an important role in a low-Reynolds-number environment the cell arrangement BHL has been used together with a prescribed addressing sequence of two neighboring actuators, V and H, forming the shape of an “L”. The two actuators are addressed in the sequence VHVHVH, which would result in a four phase cyclic motion if the actuators were to move unhinderedly. This scheme, originally proposed by one of the authors in an essay in a popular scientific journal Ricka (2010), was motivated by the “two-hinged low-Reynolds-number-swimmer”, which has been introduced in Purcell (1977) and is illustrated in Fig.5. The prescribed addressing sequence principally represents a special locally deterministic update scheme, as it considers the actual state of two locally coupled cells to determine which actuator will be addressed. The selection of the “L” is based on the update schemes introduced above (RAU1, RAU2, SRAU1, SRAU2 and DAU). Consequently, we actually mixed a local update scheme (prescribing the four-phase sequence of an “L”) with different global update schemes (selecting which “L” to address). Settings using the coupling of two cells forming an “L”, will be accounted for a cell arrangement, referred to “BHL+L”.

The different update schemes have been applied in order to investigate the role of potential intercellular signaling mechanisms on the airway epithelium for the dynamical behavior of locally interacting ciliated cells.

II.6 Local Mucociliary Interactions

Hydrodynamical interactions between adjacent ciliated cells are considered in a simplified fashion and implemented in terms of logical local decision rules induced by mucus droplets randomly seeded on the empty sites in the network. The system’s evolution is achieved by attempting to move the individual actuators. As interactions between actuators occur only if mucus is located on the activated actuator, an activated actuator can switch its state unhinderedly as long as there is no mucus hindering its oscillation. In the presence of mucus, the actuator can shift and squeeze single droplets, or stacks of droplets, depending on the available energy.

There are two possibilities that an activated actuator remains in its actual state. Either the active actuator has not enough energy to squeeze the mucus lumps on adjacent fields, or the actuator is situated in a locked configuration (see Fig.6). Note that there are always two possible locked configurations for each actuator, which we shall refer to in the following as $l_{0}$ and $l_{1}$.

If the activated actuator is not situated in a locked configuration it either flips its state by squeezing the concerned mucus on adjacent fields, or it stagnates and remains in its current state (see Fig.7). The squeezing of mucus on adjacent fields is associated with a certain change of free energy. Prior to its attempt to move, an actuator gets a certain amount of energy ascribed, which has been inverse sampled. If the actuator’s ascribed energy exceeds the required moving energy, the actuator is able to move. This way, the actuators’ energy distribution and the situation-specific change in free energy determines the probability for an actuator to flip its state by squeezing the mucus on adjacent actuators.

For computational purposes it is convenient to consider the mucus-free case ($m_{ij}^{t}=0$) as a special case, in which an actuator is always able to switch its state (see Fig.8).

Intercellular coordination, or rather the emergence of global spatio-temporal patterns, are caused by stagnating actuators as well as locked configurations, as in these situations the activated actuator has to adjust its state according to the local state and mucus distribution.

II.7 Boolean Network Representation

In this study the mucociliary dynamics on the epithelium are represented in terms of an adaptive boolean network. Ciliated cells are imported in terms of boolean actuators and represent the nodes of the network. The links between nodes represent mucociliary interactions, which we formulate in terms of boolean functions.

Formally, a boolean network consists of $N$ elements $\{\psi_{1},\psi_{2},...,\psi_{N}\}$, each of which is a binary variable $\psi_{i}\in\{0,1\}$ representing a node in the network. In general Aldana et al. (2003), the value of a node $\psi_{i}$ at time $t+1$ is given as a function $f_{i}$ of its $K_{i}$ controlling elements at time $t$:

[TABLE]

The boolean function $f_{i}$ as well as the number of controlling elements $K_{i}$ may be different for each node. The dependence on node $i$ explicitly denotes that the set of controlling nodes with indices $\{j_{1},j_{2},...,K_{i}\}$ generally varies from one node to the other.

In our case, $\psi_{ij}$ at time $t+1$ is given as a function $f_{ij}$ of the local state configuration $\{\psi_{pq}\}_{pq\in n_{ij}}$ and the local distribution of mucus droplets $\{m_{pq}\}_{pq\in n_{ij}}$ at time $t$:

[TABLE]

The neighborhood $n_{ij}$ of the actuator $ij$ is illustrated in Fig.2. The function $f_{ij}$, i.e. representing the local mucociliary interactions, imports the actuators’ available energy and will be specified in Sec.II.8.

The oscillatory driving force of ciliated cells results in the permanent attempt of a node for state reversal and is represented by the logical XOR-function:

[TABLE]

Here $c_{ij}^{t}\in\{0,1\}$ represents the local couplings among the nodes. If $c_{ij}^{t}=1$, that is if the actuator is capable of displacement, then the XOR-function inverts an actuator’s actual state $\psi_{ij}^{t}$. On the other hand, if the actuator is situated in a locked configuration (Fig.6) or stagnates due to an excessive mucus load (Fig.7), then $c_{ij}^{t}$ gets false/zero and consequently, remains in its current state.

The actuator is capable of displacement if it’s not blocked and if the free energy $u$ available to the actuator is larger than the work $\Delta w$ needed to displace the mucus. Thus, we write the condition for a state alternation as:

[TABLE]

$l_{0}$ and $l_{1}$ denote the two locked state configurations (Fig.6) and $\Theta(x)$ the Heaviside step function, such that $\Theta(0)=1$.

In order to clarify the state alternation condition we shall reconstruct the intuitively formulated local interactions, which have been illustrated in Fig.6-Fig.8:

• •

Fig.6: In locked configuration case, either $l_{0}=1$ or $l_{1}=1$ and consequently, NOR$(l_{0},l_{1})=0$. As $c_{ij}^{t}\sim\text{NOR}(l_{0},l_{1})$, the activated actuator remains in its current state.

• •

Fig.7: In hindered motion case NOR$(l_{0},l_{1})=1$ and thus the condition for state alternation reads (according to Eq.4): $u\geq\Delta w~{}$.

• •

Fig.8: Obviously, mucus-free actuators are not locked and consequently: $l_{0}=l_{1}=0$ $\implies NOR(l_{0},l_{1})=1$. As there is no mucus to redistribute the required work vanishes ($\Delta w=0$) and thus, the state alternation occurs spontaneously.

It is important to realize that the distribution of mucus droplets determines the topology of the network, as the update of mucus-free nodes is not affected by their local environment. As soon as an actuator gets occupied by a mucus droplet, it gets functionally connected to its neighboring nodes, as its subsequent state depends on the local state and mucus configuration. Consequently, in our network not only the network’s state is exposed to dynamics, but also the network’s topology, as the (re-)distribution of the droplets depends on the state transitions of the nodes and vice versa. Networks exhibiting such a feedback loop between the network’s state and its topology are called coevolutionary or adaptive networks Gross and Blasius (2008).

II.8 Actuator-Energy and Mucus Relaxation

The interactions are quantified using a simple scheme intended to model in a crude fashion the transient entropic elasticity and relaxation of entangled mucin chains. We view the actuator as a piston acting against the pressure $p$ exerted by mucus droplets contained in a certain volume $V$. For simplicity we assume that the interactions are mediated only through the set of $n$ direct neighbors of the targeted site. As a result the volume to be compressed by the action of the piston is $V=n\Delta V$ (see Fig.9).

By moving the actuator, this volume changes by $-\Delta V$ to $V^{\prime}=(n-1)\Delta V$, what requires the work $\Delta w=p\Delta V$. (Strictly speaking we should consider the pressure difference $\delta p=p_{f}-p_{b}$ between the front and back of the actuator. We apologize for forgetting the back pressure $p_{b}$ in the present set of simulation data, what fortunately does not change the main conclusions.)

To determine the pressure $p$ we employ the standard thermodynamic relation: $p/T=\partial S(V)/\partial V$, where $S(V)=k_{B}\text{ln}(W_{V})$ is the Boltzmann entropy of the mucus droplets enclosed within $V$. Replacing the differentials by discrete differences $\partial S\rightarrow\Delta S=k_{B}\text{ln}(W_{V^{\prime}}/W_{V})$ and $\partial V\rightarrow-\Delta V$ yields

[TABLE]

Subsequently, we set $k_{B}T=1$, i.e., energy is measured in units of $k_{B}T$. A plausible expression for the multiplicity (”thermodynamic probability”) $W_{n}$ can be deduced as follows: Random deposition of $M$ mucus droplets on $n$ available sites is equivalent to rolling an $n$-sided dice. Thus, the result of the deposition of $M$ droplets is a sample from the multinomial distribution, with equal probabilities for hitting a field unoccupied by an actuator. Thus, the multinomial coefficient $C_{m}=M!/(m_{1}!m_{2}!\cdot\cdot\cdot m_{n}!)$ is a good candidate for the multiplicity in Boltzmann entropy. ($m_{i}$ denotes the number of droplets deposited on a site $i$.) However, after the deposition, prior to an attempted move of an actuator, we allow the distribution of droplets to relax to a ”thermodynamic equilibrium”, i.e., into a state of maximum multiplicity, where the droplets are most uniformly distributed on the available sites, so that $|m_{i}-m_{j}|\leq 1$, for $i\neq j$. (In other words, we invoke the Maximum Entropy Principle Hanel et al. (2014).) Therefore, we define the multiplicity involved in Eq. 5 as W_{n}=\text{min}\big{(}m_{1}!m_{2}!\cdot\cdot\cdot m_{n}!\big{)} and W_{n-1}=\text{min}\big{(}m_{1}^{\prime}!m_{2}^{\prime}!\cdot\cdot\cdot m_{n}^{\prime}!\big{)}, where $m_{i}$ and $m_{i}^{\prime}$ are subject to the constraints $\sum_{i=1}^{n}m_{i}=\sum_{i=1}^{n-1}m_{i}^{\prime}=M$.

To complete the specification of the mucociliary interactions we must specify the actuators energy $u$. At this point we introduce an additional stochastic element, assuming a certain distribution $f_{U}(u)$ of the actuator energy. Prior each attempt an actuator’s energy $u$ is obtained by reverse sampling, according to: $u=F_{U}^{-1}(r)$, where $r$ represents a uniformly distributed random number and $F_{U}$ the cumulative energy distribution. Since a cumulative distribution is a monotonous function of its argument, the sampling can be concisely included into the boolean update equation for $c_{ij}^{t}$ as:

[TABLE]

In the present simulations we used (for certain “historical reasons”) the simple expression

[TABLE]

where $u\geq 0$ and $\epsilon$ parameterizes the actuators’ energy supply. This function is illustrated in Fig.10.

Finally, we summarize the update of an actuator’s state and its local mucus configuration in the form of a flow chart in Fig.11.

III Simulations and Results

III.1 Initial State

Prior to simulation the state of the virtual epithelium has to be initialized. First, an array of $N$ actuators with random states is generated, i.e. $\Psi^{0}$, where $\psi_{ij}\in\{0,1\}$ are uniformly distributed. Second, as the fraction $f$ represents the proportion of unciliated cells, $f\cdot N$ randomly chosen sites are set to $NAN$. A certain amount of mucus droplets is then seeded on a randomly chosen site (uniformly distributed). For the sake of comparability between open and toric boundary conditions the amount of mucus droplets was held constant during the simulation. Therefore, when applying open boundaries the mucus droplets leaving the modeling area were refed on a randomly chosen site. This way, we principally assumed the mucus excretion to be proportional to the mucus transport, which would of course need some kind of internal regulation mechanism in a real trachea.

III.2 Parameter Study

In order to assess under which circumstances our model self-organizes towards a self-cleaning virtual epithelium and how it reaches its properly functioning states dynamically, the influence of the model parameters, introduced in the previous chapter, on the network dynamics has been studied. The main parameter study encompasses the variation of the six model parameters in the following ranges. All cell alignments (USL, UHL, BHL, BHL+L), all update schemes (DAU, RAU1, RAU2, SRAU1, SRAU2) and all boundary conditions (OP, HC, VC, TO) have been used. The amount of mucus lumps has been varied in the range of 0.5%, 1%, 2%, 4%, …, 256% of the total number of actuators in the grid (an amount of 4% corresponds to $0.04\cdot N$ mucus droplets). The amount of unciliated cells has been set to 0%, 5%, 10%, 15%, 20%, 25% and finally, the energy parameter has been set to 0, 0.25, 0.5, 0.75 and 1. Consequently, our main parameter study encompasses 24’000 simulation runs, which have been iterated for $10^{5}$ time steps using a grid size of 50$\times$50 cells. These simulations can be seen as a starting point of further simulations we conducted. Each simulation run is characterized by its corresponding parameter setting. A specific parameter set is denoted as a combination of the form: (grid size, cell alignment, update scheme, boundary condition, mucus amount, energy parameter, amount of unciliated cells).

III.3 Observables

To provide a qualitative impression of the self-organizing character of the network model and to illustrate the meaning of the chosen observables, we present first a typical simulation run. The parameters have been set to (50$\times$50, USL, DAU, TO, 16%, 0.5, 0%). Fig.12 shows the state of the network at three different stages of the self-organization process. Fig.12 represents a substitute for the temporal evolution of the network’s state, which is best visualized by movie_S1 (see supplemental material sup ). Actuators are colored in dark gray if $\psi_{ij}^{t}=1$ or in bright gray if $\psi_{ij}^{t}=0$. Mucus lumps are shown in white. Fig.12a shows the initial network state displaying the randomly generated initial configuration of states and mucus lumps. Fig.12b visualizes the state of the network after 100 iterations, representing an intermediate stage of the self-organization process, as the emergence of global order becomes clearly visible. Finally, Fig.12c shows the network state after 1500 iterations. At this stage the fascinating self-organization process has almost completed, and the actuators finally behave strongly coordinated at a global scale.

By examining the successive network states in movie_S1, one can observe that the particles get initially moved around disorderly. Quickly, the actuators start to cooperate by adjusting their oscillations to the oscillations of the surrounding actuators until the particles get efficiently transported into a well defined direction - what we call self-organized transport. Consequently, in our models the self-organization process is actually twofold. As in accordance with the emergence of spatio-temporal patterns, the virtual epithelium exhibits self-organized transport. This co-evolution of the network state and its associated mucus transport has been quantified in terms of several observables, which are introduced in the following.

III.3.1 Mucus Transport Velocity

For each actuator at the position $ij$ at time $t$ we assign a local mucus velocity in terms of the local displacement of the center of mass (CM), which we denote as $\vec{v}_{ij}(t)$ and has been calculated (if $m_{ij}^{t}>0$) according to:

[TABLE]

$\vec{v}_{ij}(t)$ measures the redistribution of mucus droplets in the neighborhood $n_{ij}$ of the activated actuator at $ij$ in terms of the local displacement of the CM. $\vec{r}_{ij}-\vec{r}_{pq}$ denotes the distance vector pointing from the activated actuator at ${ij}$ to the neighboring actuators at $pq$. This relative distance vector gets weighted by the redistributed mucus droplets $(m_{pq}^{t+1}-m_{pq}^{t})$.

In order to quantify the global transport velocity $\vec{v}_{g}(t^{\prime})$ (remind that $t^{\prime}\doteq t/N$) we calculate the mucus-weighted average over the whole array of actuators of the formerly defined local velocity of the CM in Eq.(8) according to:

[TABLE]

$\vec{v}_{g}(t^{\prime})$ represents an approximation to the actual average mucus transport velocity, as droplets moving more than one field in one single timestep are neglected. As these movements only happen due to fluctuations originating from mucus relaxation, $\vec{v}_{g}(t^{\prime})$ represents a good measure for the area averaged mucus transport velocity.

Fig.13 shows the temporal evolution of the global mucus transport speed $|\vec{v}_{g}(t^{\prime})|$ for an ensemble consisting of 100 simulation runs differing only by their initial state. Each ensemble member corresponds to a simulation run for which the parameter setting (50$\times$50, BHL+L, RAU1, OP, 10%, 1, 0%) has been used.

The curve shows the temporal evolution of the ensemble mean (solid line) and its standard deviation (shaded area) and displays the typical saturation-like behavior of the average mucus transport speed, which has been observed for each parameter setting showing a self-organizing behavior. Accordingly, we can define the initial mucus transport velocity $\vec{v}_{g0}\doteq\lim\limits_{t^{\prime}\to 0}\vec{v}_{g}(t^{\prime})$ and the final mucus transport velocity $\vec{v}_{g\infty}\doteq\lim\limits_{t^{\prime}\to\infty}\vec{v}_{g}(t^{\prime})$. Consequently, the global average transport velocity can be expressed as

[TABLE]

where $\Delta\vec{v}_{g}(t^{\prime})$ can be seen as the effectively self-organized mucus transport velocity.

III.3.2 Mucus Transport Direction

We observed that some parameter sets drive the model towards well organized attracting network states exhibiting non-negligible area averaged transport speeds. However, some of these well organized states exhibit unrealistic velocity fields.

The network state shown in Fig.14

(upper panel) and its corresponding (temporally averaged) velocity field (lower panel) has been generated by applying the parameter set (50$\times$50, UHL, DAU, OP, 64%, 1, 0%). It can be seen that even if the network state is well organized the corresponding velocity field appears to be rather disordered. A closer look reveals that the network state as well as its corresponding velocity field is divided into two parts. The right half transports the mucus with a tendency to the right, while the transport on the left half tends to the left. The global average velocity amounts to $\vec{v}_{g\infty}=(-0.14,-0.03)$ cells/it.

In order to classify parameter sets generating such odd velocity fields as malfuctioning, we measured the spread of the (temporally averaged) local velocity fields in terms of $\langle\cos\theta\rangle\doteq\langle\cos\theta_{ij}\rangle_{ij}$, where $\langle...\rangle_{ij}$ indicates spatial averaging. We defined $\cos\theta_{ij}$ as:

[TABLE]

where $\langle...\rangle_{t^{\prime}}$ indicates temporal averaging (over the last $10^{3}$ iterations of each simulation). For the velocity field shown in Fig.14 $\langle\cos\theta\rangle$ amounts to 0.18 indicating not properly directed transport.

III.3.3 Transient Time

Since $\Psi^{t^{\prime}}$ represents the state of the Boolean network at time $t^{\prime}$ containing $N$ Boolean variables, the set of all 2N possible network states forms the state space of the Boolean network. The successive network states $\Psi^{0},\Psi^{1},\Psi^{2}...$ form a trajectory in the state space. In the case of deterministic Boolean networks, a network sooner or later reaches a state, which has been reached before (due to the finite state space) and consequently, enters a cycle consisting of a subset of states of the state space, which is called an attractor. A more general definition for dynamical systems says that an attractor is a set of states to which the system evolves after a long enough time Greil (2009). The transient time $\tau$ is the number of states a network undergoes (starting from an initial state), before it reaches an attractor Wuensche (1998); Gershenson (2004b); Greil (2012).

Fig.15 depicts the temporal evolution of two ensembles consisting of 100 ensemble members differing only by their initial condition. The brighter band corresponds to the curve presented in Fig.13 showing the evolution of the average mucus transport speed being generated with the parameter setting (50$\times$50, BHL+L, RAU1, OP, 10%, 1, 0%). The darker curve has been generated with exactly the same parameter settings apart from the update scheme. Instead of RAU1 the update DAU was used. Both ensembles show the typical saturation-like temporal evolution of the average transport speeds, which reflects the capturing of the system dynamics by attractors. We calculated the transient time based on the transient behavior of the average transport speed, which roughly follows an exponential behavior: $|\vec{v}_{g}(t^{\prime})|=|\vec{v}_{g0}|+\exp(-t^{\prime}/\tau^{\prime})\cdot(|\vec{v_{g\infty}}|-|\vec{v_{g0}}|)$ . The transient time $\tau$ has been defined as $\tau\doteq 3\cdot\tau^{\prime}$. As one can notice in Fig.15 the transient time for the ensemble simulation using the deterministic update scheme DAU is roughly four times shorter than the one for which the update scheme RAU1 was used. This behavior may be more fundamental and shall be discussed later.

III.3.4 Autocorrelation

In order to characterize the coordination amongst the actuators in a particular network state, we use the spatial autocorrelation function C$(\Delta i,\Delta j,t^{\prime})=\sum\limits_{i,j}\Psi(i,j,t^{\prime})\Psi(i+\Delta i,j+\Delta j,t^{\prime})$ , where $\Delta i$ and $\Delta j$ denote the shifts into the $i$- and $j$-direction, respectively. As we would like to compare different simulation runs with respect to the emergent spatial order, we used the autocorrelogram to determine the spatial autocorrelation-length $\rho_{c}$, representing a measure for the degree of order in a network state. As most correlograms appear to be strongly elongated as illustrated in Fig.16, we determined the autocorrelation-length along the direction of maximum correlation, which is indicated by the gray plane in Fig.16. Note that in roughly 95% of all properly self-cleaning attracting states the direction of maximum correlation coincides with the direction of mucus transport.

As most correlograms show an exponential-like decrease along the direction of maximum correlation, the autocorrelation-length $\rho_{c}$ has been determined according to:

[TABLE]

Eq.12 delivers the autocorrelation-length $\rho_{c}$ along the line of maximum correlation, which is represented by the parameterization $\gamma(v)$.

III.4 Co-evolution of Spatio-Temporal Patterns and

Transport

The co-evolution of self-organized transport and spatio-temporal patterns is illustrated in Fig.17. The brighter curve shows the area-averaged transport speed $|\vec{v_{g}}(t^{\prime})|$ (according to Eq.9) illustrating the build up of self-organized transport. The darker curve represents the auto-correlation length of the network’s state at time $t^{\prime}$. As at the beginning of the simulation the mucus droplets get displaced almost erratically, the corresponding area-averaged transport speed almost vanishes. As time passes the average particle speed grows until the coordination of the actuators reaches a maximum. The three vertical dashed lines indicate the moments at which the three network states shown in Fig.12 have been recorded.

Interestingly, all model settings leading to self-organized transport have shown a similar growth behavior with respect to the global average transport speed as the one depicted in Fig.17. This typical transient dynamical behavior reflects the restricted growth of the cooperation among the actuators. This sigmoid logistic-like dynamical behavior can be observed in other studies investigating self-organizing processes as well, and appears to be omnipresent for self-organizing processes. An example is provided by Fig.2 in Niedermayer et al. (2008), in which the self-organized synchronization and wave formation in one-dimensional cilia arrays has been studied.

Within the meaning of the state space concept, the curves in Fig.17 can be seen as the network’s transient behavior. The saturation-like behavior reflects the capture of the network dynamics by an attractor. Accordingly, the strongly ordered network state shown in Fig.12 represents an attracting state.

Our model typically exhibits a co-evolution of the network’s state and its associated mucus transport. As indicated earlier, we suggest to see this two-fold self-organization in the context of adaptive boolean networks. Since mucus droplets functionally connect an actuator to its surrounding and disconnect an actuator when squeezed away, the network topology is exposed to dynamics and interrelated to the dynamics of the network’s state. Consequently, we observe not only dynamics on the network, represented by state transitions, but also dynamics of the network, caused by the transportation of the mucus droplets. This means that changes in the network’s state are affected by the network’s topology and vice versa, forming a feedback loop between the topology and the state of the network. The resulting co-evolution of the network state and the network topology represents the characteristic property of coevolutionary or adaptive (boolean) networks Gross and Blasius (2008); Rohlf and Bornholdt (2009).

III.5 Characterisitcs of Attractor States

III.5.1 Transport of Attracting States

Fig.18 presents an overview of the terminal mucus transport ($m\cdot\vec{v}_{g\infty}$), the terminal mucus transport velocity ($\vec{v}_{g\infty}$) and the initial mucus transport velocity ($\vec{v}_{g0}$) for each cell alignment (column-wise). The first row indicates to which cell alignment each column corresponds to and illustrates the sequences after which the actuators have been activated when the deterministic update scheme DAU has been applied. Apparently, the globally averaged terminal transport velocities have a clear preference considering their direction for each cell alignment. Generally, the preferred transport direction seems to be oppositely-oriented to the direction of the deterministic update signal. However, for a small set of the settings using BHL or BHL+L the transport shares the direction of the update signal.

III.5.2 Classification Attempt

Fig.19 illustrates how the transport capability (upper colormap) and the transport speed (lower colormap) of attractor states is related to the correlation length. The gray scale indicates the frequency density.

Considering the attractors’ function and structure a rough classification into the following four classes (C1-C4) seems appropriate: C1: nontransporting disorganized states, C2: transporting disorganized states, C3: nontransporting structured states and C4: transporting structured states. We measured the transport capability in terms of the transport speed $|\vec{v_{g\infty}}|$ and the degree of the organization of the expression patterns in terms of $\rho_{c}$. States/settings exhibiting a higher $|\vec{v}_{g\infty}|$ than 0.1 (cells/it.) are considered as transporting, while states/settings exhibting a $|\vec{v}_{g\infty}|$ of less than 0.01 (cells/it.) have been considered as nontransporting. States satisfying $\rho_{c}>2$ (cells) have been classified as organized and states with a $\rho_{c}$ of less than 1.5 (cells) as disorganized ($\rho_{c}$ of the randomly generated initial states were shorter than 1.6 cells).

An illustration of how the observables and the associated parameter values are distributed in each class, as well as a listing of the most striking conspicuities, can be found in the supplemental material sup . Here, we just note that attractor states of class 4 still contain many states exhibiting odd velocity fields, like bisected or disordered ones, which still reach considerable terminal transport speeds. Consequently, only a subset of class 4 is considered as “properly self-cleaning” attractor states, whose velocity fields shall fulfill further conditions specified in the next section.

III.5.3 Effectively Self-Organized Properly Self-Cleaning

Attractors

In order to find the parameter values allowing the network to self-organize towards properly self-cleaning attractor states, we finally classified the settings as “functioning” or ”malfunctioning”. A certain parameter set is classified as functioning, if the following four conditions are simultaneously fulfilled after an integration time of $t^{\prime}=10^{5}$.

1. The globally averaged transport speed $|\vec{v}_{g\infty}|$ is faster than 0.1 cells/it., which would roughly correspond to 20$\mu\text{m}/s$ in the real system (assuming a ciliary beat frequency of 10 Hz and the diameter of a ciliated cell being 10$\mu\text{m}$).
2. We require the velocity field being sufficiently directed by: $\langle\cos\theta\rangle>0.65$.
3. The velocity field must be self-organized and not simply being imposed by the choice of a parameter set leading to an initial tendency of the transport direction. Consequently, by requiring $|\Delta\vec{v}|\doteq|\vec{v_{g\infty}}-\vec{v_{g0}}|>0.1$, we demand that the velocity field changes (at least its direction) in the course of a simulation.
4. The auto-correlation length $\rho_{c}$ has to be longer than two cells. The number of iterations may be seen as a fifth condition considering the efficiency of the self-organization process. As thus, the transient time $\tau$ has to be shorter than $10^{5}$ iterations. The classification of parameter sets as functioning and malfunctioning was applied to the parameter study encompassing 24’000 parameter sets. 564 out of 24’000 settings ($\widehat{=}$ 2.4$\%$) have been classified as functioning. Table 1 shows how these functioning settings spread across the different cell alignments and update schemes (values are listed in $[\%]$).

Table 1 summarizes the following observations. The cell alignment USL can cope with the largest set of parameter settings and makes up $84.2\%$ of all functioning settings. $98.2\%$ of all functioning settings are made up by settings using either the cell alignment USL or BHL+L. The deterministic update scheme DAU can cope with each cell alignment and $89.6\%$ of all functioning settings were generated by this update scheme. None of the settings using the completely random update RAU2 evolves the network to a self-organized self-cleaning state. The cell alignment BHL+L (lowest row in table 1) can cope with all update schemes except the fully random update scheme RAU2. Apparently, a certain degree of local determinism is advantageous for an efficient self-organization towards properly transporting attractor states, since all functioning parameter settings either involve DAU or BHL+L.

III.5.4 Mean Wave Propagation Direction in Stereotypical

Functioning States

The set of functioning parameter values drives the model towards families, or stereotypes, of attractor states comprising very similar attractor states. In order to determine the mean direction of wave propagation in these sterotypical attractors, we employed the concepts outlined in Ryser et al. (2007) and examined sequences of sequential attracting states of one specific representative parameter set for each stereotype. A more detailed discussion of these representative stereotypical attractor states can be found in the supplemental material sup . Here, we would like to report the summarizing comparison between the direction of the mean wave propagation and the direction of transport, what can be found in Table LABEL:tab:tab2. Note that mean $k$-vectors visualized by double arrows indicate that the mean direction of wave propagation could not be determined unambiguously, which comes from the boolean nature of the model (undersampling).

III.6 Dynamical Characteristics

III.6.1 Imposed Asymmetry as a Starting Assistance

For Self-Organized Self-Clearing

In this section we discuss the potential role of imposed asymmetries as a control parameter. As shown in Table 1 functioning settings either include DAU or BHL+L, or both - all other settings don’t result in self-organized self-clearing states. Consequently, these settings may reveal a common characteristic which ultimately leads to self-organized self-cleaning. We realized that properly functioning parameter sets consistently impose an initial tendency of the transport direction. The initial area averaged transport speed $|\vec{v}_{g0}|$ can be used in order to quantify the anisotropy of the transport at time $t^{\prime}$=0. In order to estimate the expected initial transport speed at $t^{\prime}$=0 of a certain parameter setting, we used a grid consisting of 1500$\times$1500 actuators, which have been updated for a single timestep ($\Psi^{0}\rightarrow\Psi^{1}$). This has been done for an ensemble of 100 simulation runs. The ensemble mean provides an estimate for the imposed initial transport speed $|\vec{v}_{g0}|$. Fig.20 illustrates the frequency distribution of functioning (dark gray bars) as well as malfunctioning (bright gray bars) parameter sets consindering their corresponing $|\vec{v}_{g0}|$-classes (10 logarithmically spaced $|\vec{v}_{g0}|$-classes between $10^{-5}$ and 1). Note that the scale of the two ordinates differ by one order of magnitude.

It can be particularly seen that the relative frequency of functioning parameter sets increases with increasing $|\vec{v}_{g0}|$-values. Moreover, all functioning parameter sets correspond to $|\vec{v}_{g0}|$-values being faster than 0.01 cells/it.

Correspondingly, the measures $|\Delta\vec{v}|$, $\rho_{c}$, $\langle\cos\theta\rangle$ and $|\vec{v_{\infty}}|$ classifying the parameter sets tend to increase with increasing values of $|\vec{v}_{g0}|$, which is illustrated in Fig.21.

Consequently, Fig.20 and Fig.21 indicate that $|\vec{v}_{g0}|$ may trigger the model towards self-organized self-clearing.

The boxplots shown in Fig.22 (Whiskers are set in order to indicate the range covered by 95% of all values) display how the $|\vec{v}_{g0}|$-values distribute over the different update schemes and cell alignments. One can see that the highest values for $|\vec{v}_{g0}|$ are reached by settings using DAU or BHL+L. This clearly suggests that update schemes with local determinism impose an asymmetry among the local interactions, what triggers efficient self-organization towards self-clearing states.

III.6.2 Effect of The Update Scheme on The Network Dynamics

As settings using the locally prescribed four-phase sequence (BHL+L) can cope with each update scheme, we compared the dynamic behavior produced by the different update schemes by simulating an ensemble of 100 simulation runs using the parameter settings: (50$\times$50,BHL+L, $\ast$, OP, 10%, 1, 0%), where $\ast\in$ $\{$DAU, RAU1, SRAU1, SRAU2$\}$. Fig.23 shows the transient behavior for the average mucus transport speed of the four ensembles corresponding to the different update schemes. Each setting produces a saturation-like behavior reflecting the capturing of the dynamics by attractors. As one can see, settings using SRAU2 and DAU reach their attractors considerably faster than settings using SRAU1 and RAU1.

In order to characterize the state space corresponding to each setting, we counted the number of attractors and the attractor periods. It turned out that each attractor has a period of four, corresponding to the prescribed cyclic four-phase sequence. Furthermore, SRAU2 and RAU1 drive each of the 100 different initial conditions towards the same attracting state, which is shown in Fig.24. DAU and SRAU1 produce more realistic dynamics, as the 100 simulation runs were driven towards 90 and 95 different attracting states, respectively. For DAU and SRAU1 two exemplary attracting states are depicted in Fig.24.

Even if SRAU1 drives different initial conditions towards different attracting states, the attracting states strongly resemble each other, as indicated in Fig.24. DAU produces more different emerging structures. In order to verify this observation, we ran 5 simulation runs for each update scheme, which differed by their initial state. But this time we prescribed how much their initial state shall differ against a reference run. Namely, we switched the state for 1%, 2%, 3% and 4% of all actuators, respectively, and compared the networks’ state at each timestep to the reference run in terms of the normalized Hamming distance 111Oscillations of 0.5 in the normalized Hamming distance (caused by the XOR-function and the BHL+L setting) are suppressed by taking the minimum of four consecutive Hamming distances. (see e.g. Greil (2009)) see Fig.25). The average Hamming distance corresponding to SRAU2 and RAU1 decreases to zero, which means that all states reach the same attractor. On the other hand, SRAU1 and DAU produce different attractors, but the attractors produced by DAU show an almost maximum Hamming distance. This behavior reflects the ability of the system to find new attractors in case of a perturbation and to conform to changes, which is an important characteristic for living systems. Consequently, if an attractor which has been reached by applying DAU gets perturbed, the system conforms to changes and runs into a new attractor.

III.6.3 Open Boundaries: Marginal Nodes Guiding The Structure Emergence

First of all we consider a completely dense carpet of ciliated cells represented by a parquet of actuators. In this “densely ciliated carpet case” we observe that settings with open boundaries yield a structure emergence always starting at the same open boundary, from which it spreads over the whole network. This behavior is illustrated in Fig.26, where three snapshots of the network state at three different stages of the self-organizing process are shown. The graphs have been generated by applying open boundaries, RAU1 and BHL+L.

Further, we observed that settings leading to the “crystallization-process”, shown in Fig.26, only show this self-organizing behavior as long as the boundary, from which the structure spreads, is open. This means that similar network dynamics have been observed for open and horizontal cylindrical boundaries. On the other hand, if we chose toric or vertical cylindrical boundaries, for which the “structure-triggering” border is glued to the opposing border, the model displays a completely different dynamical behavior with much lower self-organized transport speeds and much less well ordered network states. The only exception is given by the settings using the square-lattice alignment, which shows self-organized directed transport under each boundary condition.

Movie_S2, movie_S3 and movie_S4 (provided as supplemental material sup ) illustrate the structure emergence from an open boundary. They have been generated by applying (50$\times$50, USL, DAU, OP, 10%, 0, 0%), (50$\times$50, BHL+L, DAU, OP, 15%, 1, 0%) and (50$\times$50, BHL+L, RAU1, OP, 15%, 1, 0%), respectively.

III.6.4 Leaders and Followers

Since we observed for the “densely ciliated carpet case” when using open boundaries that the emergence of highly ordered structures always sets in from the same boundary, the actuators at the boundary seem to play an important role for the self-organizing process and consequently, it appears that there exists a kind of hierarchy among the actuators.

This hierarchy among the nodes is caused by the underlying network topology and is primarily characterized by the distribution of in- and out-degrees. In Fig.27 we illustrate the underlying network topology for an array consisting of $5\times 5$ cells arranged in a square lattice assuming open boundaries. The middle panel in Fig.27 illustrates the network topology considering the possible pathways of the mucus droplets. The right panel illustrates the network topology considering the possible locked configurations. Arrows entering a node indicate which nodes may block its oscillation. Arrows leaving a node indicate for which nodes its state and mucus load is relevant to block the nodes the arrows are pointing at. The in-degree of a node is the sum of incoming arrows and correspondingly, the out-degree is the sum of outgoing arrows. The grayscale in the graphs indicate the value of the in- and out-degrees and therewith the prevalent hierarchy among the nodes.

If we set open boundaries, marginal nodes in our networks are hardly influenced by their adjacent nodes. It is especially important that nodes at the left and right boundary have no in-degree considering the locked-configurations-rule. Consequently, these nodes won’t adapt their oscillatory motions according to their neighbors’ motion and act as “leading nodes”.

Finally, we observed that some settings need a very regular topology with a prevalent hierarchy among the agents. If this “leadership” of a few nodes is abandoned, which happens when applying toric boundary conditions or importing unciliated cells, hardly any self-organization takes place.

III.6.5 Boundary-Driven Structure Emergence Hindered by Unciliated Cells

Unciliated cells can only take up and release mucus droplets and are thus, only passively involved in the state update of adjacent actuators. This means that cells adjacent to unciliated cells obtain a lower in-degree and therefore ascend the hierarchy of actuators. As we randomly distribute unciliated cells, “leading nodes” are no longer only found at the boundaries, but rather are distributed allover the network and consequently, the network topology becomes less regular.

It has been observed that the introduction of unciliated cells hinder the spreading of the structure emergence for those settings for which the structure spreads strictly from an open boundary. Consequently, it seems that some settings require a very regular topology in order to be able to spread allover the network. The panels in Fig.28 show three stages of a simulation for which we introduced 1$\%$ of unciliated cells and otherwise applied the same settings as in Fig.26. Unciliated cells are shown in black. It’s clearly visible that the structure emergence is hindered by the introduced unciliated cells, as the structure can’t spread further than to the first unciliated cells, seen from left. According to this simple observation one could imagine that the influence of the boundaries gets the more restricted the more unciliated cells are incorporated. Further simulations have confirmed this idea.

III.6.6 Modular Self-Organization

If we arrange the actuators in the square lattice (USL) and introduce a certain amount of unciliated cells, the array of actuators efficiently evolves to a self-cleaning epithelium. In this case the network topology considering the state update gets strongly changed, what becomes obvious when thinking of the case of sparsely distributed ciliated cells. One can imagine that a group of ciliated cells would be surrounded by unciliated cells forming thereby “ciliated islands” on an otherwise unciliated epithelium. These islands would be hardly interconnected amongst each others when concerning the state updates. How strongly these modules are interconnected depends on the density of ciliated cells. In order to illustrate this modular character introduced by unciliated cells, the topology of the square lattice has been applied together with open boundary conditions and the update DAU for an ensemble with 100 members differing by their initial state. The grid size has been set to $300\times 300$ cells. Fig.29 illustrates the temporal evolution of the corresponding velocity fields. The grayscales illustrate the locally resolved transient time $\tau_{ij}$ (in effect the number of network updates it takes for an actuator to reach a local transport speed, which almost amounts to the final area averaged transport speed). The upper panel in Fig.29 corresponds to simulations ran with a dense mat of cells (100% ciliated cells), while in the lower panel 10% of randomly distributed unciliated cells were introduced. The panels represent the ensemble average of the spatially resolved transient times. The darker the color the longer it took until a certain actuator synchronized its movement. The discussed effect of the “leading boundary” is clearly visible in the upper panel of Fig.29, as the gradient in brightness from the left to the right indicates the structure emergence, which sets in at the left boundary and spreads towards the right boundary. On the other hand, as soon as one introduces unciliated cells the self-organization process gets a modular character, as the structure emergence doesn’t start at a certain point from which it spreads allover the network, but spreads at several locations simultaneously, which is illustrated by the relatively homogeneous distribution of grayscale in the lower panel in Fig.29. The influence of the boundary is strongly restricted to the marginal nodes at both sides. The larger one chooses the grid size, or the more unciliated cells are introduced, the less important the “dominance” of the marginal actuators gets.

In conclusion, unciliated cells introduce a certain degree of modularity in the network topology as well as in the self-organizing process and therefore, the influence of the choice of boundary conditions can be neglected in the interior of the array, if the array size is chosen large enough and if a realistic amount of unciliated cells is considered.

Movie_S5 (provided as suplemential material sup ) shows the evolution of the network’s state of a simulation run generated with (100$\times$100, USL, DAU, OP, 20%, 0.5, 10%) exhibiting the typical modular self-organization process leading to modular expression patterns.

III.6.7 Transient Time vs. Network Size

In this section we point out that unciliated cells may strongly influence the dynamics of ciliated epithelia. As outlined in the previous section unciliated cells import topological modularity. Modularity is an important and promising concept, which is inter alia studied in terms of modular random boolean networks. In Poblanno-Balp and Gershenson (2011) it is claimed that topological modularity reduces the probability for damage spreading over the network, what promotes robustness.

In the following we show how unciliated cells affect the dynamics of our network model. As the square lattice is the only cell arrangement leading to a properly self-cleaning state, when introducing unciliated cells, all the results shown in this section refer to the square lattice arrangement.

The transient time has been determined for different array sizes in the range between $50\times 50$ cells2 and $1000\times 1000$ cells2 as well as for different fractions of unciliated cells (0%, 2% and 10%). The collected data points are presented in Fig.30. The lines are only guidelines for the eye. The dots correspond to the dense mat case for which all cells are ciliated. Diamonds and stars correspond to 98% and 90% ciliated cells, respectively. One can clearly see that the transient time is not only reduced when introducing unciliated cells, but varies substantially different with increasing array size. While the transient time continuously grows for an array representing a totally ciliated mat with increasing array size, the transient time runs into saturation with increasing array size if unciliated cells are present. As discussed in the former section unciliated cells import not only topological modularity, but also cause the self-organization process to be modular, as the emergence of structure and transport evolves simultaneously in different modules. This finding most probably explains the level off of the transient time with increasing array size if unciliated cells are considered. As the self-organization process takes place in a decentralized manner, it does not depend upon the size of the actuator-array. On the other hand, the transient time continues to increase with increasing array size for an array purely consisting of ciliated cells, which is related to the nature of the structure emergence for these settings. According to the upper panel in Fig.29 the self-organization process mainly starts at the left boundary, from where it spreads to the right allover the network. Accordingly, one would expect that the growth of the transient time continues with increasing network size, as it is the case in Fig.30. (The transient time primarily depends on the length of the simulated actuator-array, as the ordered structures expand from the left to the right).

As the transient time can be seen as a measure for robustness, our results suggest that the imported modularity promotes robustness.

Consequently, it seems that the modular topology imports modularity into the self-organization process, which means that the organization of a huge network gets decomposed into a simultaneous organization of submodules. These submodules can be recognized by examining the spatial structure of the network’s state, as it has been observed that the modularity of the expression patterns clearly depend on the density of ciliated cells. The auto-correlation length of attracting network states has been determined for settings differing by the amount of ciliated cells as well as the boundary conditions. As we are interested in the spatial structures of the expression patterns, the grid size has been chosen to 200$\times$200 cells. Fig.31 shows the auto-correlation length $\rho_{c}$ (cells) vs. the relative amount of unciliated cells (%) for open and toric boundary conditions.

It can be clearly seen that the auto-correlation length, decreases with an increasing portion of unciliated cells. Furthermore, the auto-correlation length is roughly given by the mean distance of unciliated cells. Consequently, unciliated cells may play a further role in the self-organization process on the ciliated airway epithelium, as the appearance of the previously described patch-work character may be strongly influenced by the distribution of unciliated and ciliated cells.

IV Discussion

The aim of this study is twofold: On the one hand, we want to make the self-organized spatio-temporally coordinated ciliary beat patterns as well as the self-organized fluid transport across multiciliated epithelia plausible. We suggest that the cooperation among ciliated cells emerges from locally interacting oscillating cilia bundles belonging to different ciliated cells. As our goal was to keep our model as simple as possible, we present a virtual self-cleaning epithelium model based on symmetrically interacting two-state actuators, which we formulate in terms of an adaptive boolean network. In the framework of adaptive boolean networks the oscillatory motion of ciliated cells can be represented by “blinking nodes” and discrete mucus droplets establish the local interactions and therefore the network’s topology. In Sec.III.4 we demonstrate the coevolution of the network’s state and its topology, which is a characteristic property of adaptivity and represents the self-organized coevolution of ciliary beating patterns and associated fluid transport.

On the other hand, we report our insights to our system’s dynamics we gained by conducting parameter studies. In the following we discuss the observed effects of the update scheme, the boundary conditions and the amount of unciliated cells on the dynamics of our network model.

As we formulate our epithelium model in terms of a discretized asynchronous multi-agent network, the question of how to update the network arises. The introduced update schemes are meant to represent different possible intercellular signaling mechanisms (membrane potentials and calcium waves). Only settings using either the deterministic asynchronous update (DAU) or the alignment BHL+L, or both, guide the network towards properly self-cleaning states (Table 1). Settings using BHL+L or DAU introduce a local recurrent temporal coupling among adjacent actuators inducing an asymmetry of the local interactions. Fig.20 suggests that the prevalent asymmetry, such as the initial speed associated with a specific parameter setting, increases the probability for efficient self-organization towards self-cleaning states. In any case none of the roughly 11’000 simulations exhibiting an initial transport speed of less than $0.01$ cells/it self-organizes efficiently towards a properly self-cleaning state. The fundamental role of asymmetric interactions on extended systems has been discussed previously. Asymmetry induced effects on the synchronization process of a pair of coupled fields have been reported in Boccaletti et al. (2005), where it has been particularly argued that small changes in the asymmetry of the interactions could be used as an efficient way to synchronize or desynchronize the dynamics, as well as select the main statistical properties of the synchronized motion in ensembles of interacting units and consequently, may have relevant consequences in natural systems. In Ghorbani and Najafi (2016) the synchronization process of a ciliary chain attached to a cylindrical surface has been investigated. Each cilium is modeled in terms of a small sphere moving along an elliptic trajectory. It has been shown that an asymmetry in their orbits trigger the emergence of metachronal waves. Symmetrical settings have not shown any correlations in their beating patterns, what compares well to our results.

Furthermore, the application of the deterministic update scheme (DAU) seems to generate less well ordered attracting network states. This behavior is exemplary illustrated in Fig.24 showing a perfectly ordered attracting state for the random asynchronous update (RAU1), slightly less well ordered network states for the semi-random update schemes (SRAU1 and SRAU2) and finally, the least ordered attracting states for the deterministic update scheme (DAU). All settings using the deterministic update scheme have consistently generated patterns showing “defects” perpendicular to the direction of the update scheme and maximum correlation along the direction of the update. Fig.16 represents an autocorrelogram of an extreme case: maximum correlation is found into the direction of the primary update direction (from right to left), which is most probably caused by the local temporal coupling among actuators, while there seems to be almost no local temporal coupling from the top to the bottom, what leads to typically elongated autocorrelograms if DAU is applied.

Finally, the less strict organization the network generates when using DAU may lead to more flexible dynamics as in the case of a perturbation the system does not simply recover its original attracting state, but conforms to the changes by running into a completely different attractor (see Fig.25).

The effects of different update schemes on the dynamics of multi-agent systems are still being investigated. In Cornforth et al. (2005) six different update schemes have been applied on one-dimensional cyclic cellular automata to compare the resulting dynamics. It has been concluded that deterministic update schemes confer a degree of flexibility upon the system dynamics, what compares well to the observed conforming character of our model settings using a deterministic update scheme. Consequently, so far, the findings show evidence that in various asynchronous processes leading to self-organization, a deterministic update scheme leads to more realistic dynamics (flexibility and robustness) and may therefore be favored by evolution Gershenson (2004a). Recall, however, that “deterministic” does not mean that actuators would displace the mucus in a pre-determined direction. Rather, the transport direction evolves through the interplay between the update asymmetry and the largely stochastic interactions.

Finally, we would like to point out the possible effect of unciliated cells on the dynamics on ciliated epithelia. The topology of our network model is primarily given by the formulation of the local interaction rules, the choice of the boundary conditions and the amount of unciliated cells. As we have outlined in Sec.III.6.6 the amount of unciliated cells introduces a certain degree of topological modularity. The topological modularity in turn causes a modular self-organization, which means that the self-organization does not start at a specific point or boundary on the grid - as it has been observed for completely dense mats of ciliated cells - but starts in each module simultaneously. This modular character of the self-organization leads to the size-independence of the transient time, reported in Sec.III.6.7. Consequently, modularity may provide robustness even to networks as large as the human ciliated airway epithelium consisting of more than $10^{9}$ cells Mercer et al. (1994), as perturbations quickly fade away. Furthermore, the modular topology leads to modular expression patterns, the size of which are roughly given by the mean distance of unciliated cells (as shown in Fig.31). Finally, the finding of a modular self-organization caused by the underlying modular topology providing a highly robust patch-work amongst actuators, provides a consistent explanation of the modular expression patterns previously reported in experimental studies aimed at the quantitative description of the modulation wave fields on the tracheal epithelium Ryser et al. (2007).

As very recently pointed out in Dey et al. (2017), theoretical studies investigating the collective dynamics of hydrodynamically interacting cilia have, so far, usually considered homogeneous carpets of cilia. In Dey et al. (2017) the role of a spatial heterogeneous ciliary distribution on coherent ciliary beating using one dimensional arrays of cilia represented by rowers has been investigated. It is particularly shown that the phase coherence of random clustered distributions of rowers are less sensitive to variations of the number density than (homogeneously) random distributions of rowers. This finding might be seen as another specific dynamical phenomenon improving robustness by an underlying modular (network) topology.

We conclude that an intercellular signaling mechanism is probable on ciliated epithelia, as deterministic update schemes drive the model towards robust self-organized states, which can still conform to changes. We suggest that the patchy expression patterns of the modulation wave field observed on real ciliated epithelia may be the result of the underlying modular topology, which is primarily formed by the distribution of ciliated and unciliated cells. This patch-work character among ciliated cells may be highly robust due to a modular self-organization, which prevents perturbations to spread over the whole network. Furthermore, the boundary conditions may become irrelevant on epithelia being either large enough or having a low amount of ciliated cells.

We close this study by hypothesizing that the modular organization of the dynamics on ciliated epithelia may be seen as a robust size-independent construction plan of nature, which leads to properly self-cleaning airways in organisms being as small as new born mice as well as in adult giraffes.

V Classification Attempt

In the following we shall classify the states reached after $10^{5}$ iterations (and their corresponding parameter settings) into four classes:

• •

Class 1 (C1): nontransporting ($|\vec{v}_{g\infty}|<0.01$) disorganized state ($\rho_{c}<1.5$)

• •

Class 2 (C2): transporting ($|\vec{v}_{g\infty}|>0.1$) disorganized state ($\rho_{c}<1.5$)

• •

Class 3 (C3): nontransporting ($|\vec{v}_{g\infty}|<0.01$) structured state ($\rho_{c}>2$)

• •

Class 4 (C4): transporting ($|\vec{v}_{g\infty}|>0.1$) structured state ($\rho_{c}>2$)

Fig.S1 shows how the correlation lengths of the 24’000 randomly generated initial states are distributed and can be seen as a justification of the chosen classification into organized and disorganized attracting states.

The histograms in Fig.S2 illustrate the distribution of $\rho_{c}$-, $\langle\cos\theta\rangle$-, $m|\vec{v}_{g\infty}|$-, $|\vec{v}_{g\infty}|$- and $|\vec{v}_{g0}|$-values (row-wise) in each class (column-wise).

Accordingly, the histograms in Fig.S3 illustrate how the parameter values are distributed in each class.

In the following we list the most striking conspicuities of each class and explain some of the mechanisms accounting for the classification of a certain parameter set.

Class1 : nontransporting disorganized states

The very low $\langle\cos\theta\rangle$-values indicate chaotic velocity fields.

Settings using BHL+L are scarcely found in this class, which is due to their high initial transport speeds generated by the pre-coordinated cyclic four phase motion.

Two thirds of the settings in this class are made up by settings using either RAU1 or RAU2, which is most probably because these update schemes do not impose any asymmetries among the local interactions.

Class2 : transporting disorganized states

In some cases the states reach high $\langle\cos\theta\rangle$-values. This means that organized velocity fields are possible even if the states are disorganized.

Two thirds of the states in this class are generated by BHL+L, which is due to their high initial speeds ($|\vec{v}_{g0}|>0.1$ cells/it.).

The terminal transport speeds in this class are higher than 0.1 cells/it. and less than 0.3 cells/it.

This class is interesting because one might wonder if it is possible that transport can develop without associated structure emergence. A detailed look into the simulations showing a considerable increase of transport speed, has shown that spatial correlation always slightly increases in tandem during the course of a simulation.

Class3 : nontransporting organized states

This class is characterized by very small terminal transport speeds as well as disorganized velocity fields (low $\langle\cos\theta\rangle$-values), but nevertheless highly organized states (two states even reach a correlation length of $\approx$10 cells).

The initial transport speeds are remarkably high.

Three out of four states in this class are generated by settings using BHL+L.

Further, high mucus amounts as well as low $\epsilon$-values seem to constitute the majority of states/settings.

Visual examination of the simulations has shown that BHL+L settings can generate ordered states before the actuators’ motion gets completely blocked by an overload of mucus. This usually can happen with moderate mucus amounts and $\epsilon<1$, or less frequently, with high mucus amounts when $\epsilon$ is set to 1.

Settings using USL and $\epsilon=1$ can generate structures. However, during the course of a simulation the mucus particles accumulate until reaching an unnatural state, which shifts piles of mucus back and forth resulting in no net movement.

Settings using UHL reach fully blocked states as well as structured states exhibiting velocity fields which are divided into two parts (similar to the example shown in Fig.14 in our manuscript) and whose globally averaged transport speeds vanish.

Class4 : transporting organized states

In this class we observe the highest values for each observable. The majority of parameter sets either involves USL or BHL+L together with the deterministic update DAU.

Further analysis of this class has shown that the direction of maximum correlation of roughly 90% of the attractor states coincides with the transport direction.

This class still contains many states exhibiting strange velocity fields, which either appear disordered or bisected. Especially bisected states can nevertheless reach considerably high globally averaged transport speeds, while their velocity fields are characterized by low $\langle\cos\theta\rangle$-values.

VI Superior States/Settings

In Fig.S4 we present the state, its corresponding autocorrelogram and its associated (temporally and spatially averaged) velocity field of the state exhibiting maximum transport (maximum $m\cdot|\vec{v_{g\infty}}|$) for each cell alignment (row-wise). Correspondingly, Fig.S5 and Fig.S6 show the state, its autocorrelogram and its velocity field of the settings exhibiting maximum transport speeds and displaying the highest correlation lengths, respectively. Table 1 lists the applied parameter sets and Table 2 the associated observable values.

VII Discussion of Stereotypical Self-Organized Self-Cleaning Attractors

In the following we discuss the properties of selected representatives of “families” of attractor states classified as functioning in the accompanying manuscript.

VII.1 Wave Field Analysis

In order to determine the wave propagation direction in attractor states we shall employ the concepts outlined in great detail in Ryser et al. (2007). The algorithms presented in Ryser et al. (2007) were developed for the analysis of experimentally collected videos of the mucus modulation wave fields generated by the coordinated movements of underlying tracheal cilia.

Here, we adopt the definition of the spatial power spectrum P$(k_{j},k_{i})$ as well as of the space-time correlation function C$(\Delta j,\Delta i,\Delta t)$.

The spatial power spectrum P$(\vec{k})$ represents the distribution of the wave-vectors $\vec{k}=(k_{j},k_{i})$ and is defined as:

[TABLE]

In Eq.S1 $F(k_{j},k_{i},f)$ denotes the distribution of complex Fourier amplitudes, which was determined by applying a three-dimensional fast Fourier transform (FFT) to a sequence of attractor states ($\Psi(i,j,t^{\prime})$, where t’ $>\tau$).

The space-time correlation function C$(\Delta j,\Delta i,\Delta t)$ is given as:

[TABLE]

where $t^{\prime}>\tau$. The velocity of the “mean wave” $\vec{v}_{wave}$ can then be measured by

[TABLE]

where $\vec{\sigma}=(\sigma_{j},\sigma_{i})$ denotes the displacement of the correlation maximum with respect to the origin $\Delta j=0,~{}\Delta i=0$.

VII.2 USL

VII.2.1 Stereotype 1

The parameter setting (50$\times$50, USL, DAU, HC, 256%, 0.25, 0%) generates the attractor state with the highest transport rate, which amounts to $m\cdot|\vec{v}_{g\infty}|\approx 4\cdot 10^{3}$ (#$m\cdot$cells/it.). The terminal velocity reached amounts to: $\vec{v}_{g\infty}=(-0.69,0)$.

A snapshot of the considered attractor state is shown in Fig.S4. Once reached this pattern remains stable. Strictly speaking, assuming the state shown in Fig.S4 is represented by $\Psi(i,j)$, the reached attractor state alternates bettween $\Psi(i,j)$ and its (logical) negation $\overline{\Psi(i,j)}$. Accordingly, once rached the terminal transport velocity remains constant, as is illustrated in Fig.S7.

Consequently, a sequence of attractor states appears like a ’binary wave’ propagating horizontally. However, this binary wave may actually represent a standing wave. This leads to a temporally constant spatio-temporal correlation function C($\Delta i,\Delta j,\Delta t$) and thus consisits of a sequence of the same autocorrelogram, which is shown in Fig.S4.

Consequently, this attractor state may either represent a (binary) standing wave or it might be constituted of waves propagating either to the right or to the left, which can not be determined unambiguously.

Finally, all functioning attractor states reached by USL, for which either an open boundary or horizontal cylindric boundary conditions have been used, closely resemble each other and particularly reach a stable state representing either a horizontally propagating or a standing wave.

VII.2.2 Stereotype 2

The parameter setting (50$\times$50, USL, DAU, TO, 64%, 0.25, 0%) generates another representative of a stereotypical attractor state exhibiting well-directed efficient transport. The pattern reached by this setting does not remain constant over time. Instead, it constantly reshapes while keeping its modular appearance. A snapshot of a particular pattern of this complex attractor is provided by Fig.S8a . Fig.S8b shows the corresponding mean auto-correlogram C$(\Delta i,\Delta j,0)$. Similar to the first stereotypical pattern, the autocorrelogram is again elongated in parallel to the transport, which is directed from right to left: $\vec{v}_{g\infty}=(-0.5,0)$. However, this stereotype exhibits a spatially as well as temporally decreasing correlation function C$(\Delta i,\Delta j,\Delta t)$. The temporal decrease of the correlation function is illustrated in Fig.S8c, which shows the profile of the correlation function along the horizontal (given by $\Delta i=0$) for increasing time lags. Before the correlation is lost the maximum of the correlation function shifts its position along the horizontal to the right, as indicated by the arrow. We observed that a higher amount of “interaction particles” (or a lower value for $\epsilon$) causes the mean wave to move faster to the left as well as a faster de-correlation.

Finally, Fig.S8d shows the transient behaviour of the average transport speed $|\vec{v}_{g}|(t)$. It can be seen that $|\vec{v}_{g}|(t)$ keeps varying around a typical average transport speed, which reflects the permanent re-organization of its associated state.

VII.3 UHL

VII.3.1 Stereotype 1

A representative of the first stereotype when using the cell alignment UHL gets generated with the parameter set (50$\times$50, UHL, DAU, VC, 128%, 1, 0%). A snapshot of the corresponding attractor state is shown in Fig.S9. Once reached this state remains stable (besides its alternation between the illustrated state and its logical negation) and its corresponding terminal transport velocity amounts to $\vec{v}_{g\infty}=(-0.26,-0.28)$. The binary nature of the considered attractor state does not allow us to determine the propagation direction of the mean wave. Consequently, the dominating wave-like binary structure either represents a wave propagating along or against the transport direction, or it might actually represent a standing wave.

VII.3.2 Stereotype 2

The second stereotype gets generated with toric boundary conditions and otherwise similar parameter values, an example is provided by: (50$\times$50, UHL, DAU, TO, 256%, 1, 0%). The attractor reached by this setting is characterized by mono-phase areas gliding over the network. These areas/gliders typically have a triangular-like shape, as indicated by the snapshot provided in Fig.S10a.

Even if the terminal transport is slow: $\vec{v}_{g\infty}=(-0.06,0.11)$, the correspoding velocity field is well-directed (see Fig.S10b).

The strong pixelation of the correlation function made it difficult to track the maximum of the spatio-temporal correlation function. Therefore, the direction of the weighted mean wave vector provides a more accurate estimate for the wave propagation direction. The weighted mean wave vector amounts to $\langle\vec{k}\rangle=(-0.05,0.05)$. The mean autocorrelation function and the relative power spectral density P$(\vec{k})$ is shown in Fig.S10c and Fig.S10d, respectively.

VII.4 BHL

Settings using the bidirectional hexagonal alignment generate very similar functioning attractor states with vertical cylindric and toric boundaries. The best transporting state classified as functioning is generated by (50$\times$50, BHL, DAU, TO, 256%, 1, 0%). A snapshot of the generated attractor is shown in Fig.S11. The state as well as sequences of states look rather chaotic, which is reflected in the broad distribution of $\vec{k}$-vectors. However, a slight tendency of waves propagating to the top right can be observed by the bare eye. Correspondingly, the expectation value amounts to $\langle\vec{k}\rangle=(0.04,0.03)$. Even if the wavefield appears rather disordered, the velocity field is well organized, since the (temporally averaged) local velocities tend to head to the right (see Fig.S11b). The terminal transport velocity is rather slow and amounts to: $\vec{v}_{g\infty}=(0.2,0.0)$.

VII.5 BHL+L

VII.5.1 Stereotype 1

The usage of open or horizontal cylindric boundaries and the update scheme RAU1 leads to a crystal-like structure. A particular example generated by the parameter set (50$\times$50, BHL+L, RAU1, OP, 256%,1,0%) is shown in Fig.S12.

Since settings including BHL+L impose a prescribed four-phase cycle on two actuators forming an “L”, we do not further have to deal with binary images and can exploit the accessibility of four phases for the analysis of the wavefield.

The particular stereotypical attractor corresponding to the snapshot shown in Fig.S12 exhibits a plane wave, whose propagation direction can not be determined unambiguously, since in this particular state the exhibited wave remains binary, despite the accessibility of four L-phases (L’s along the horizontal remain sperated by phase shifts of $\pi$). The terminal transport velocity amounts to: $\vec{v}_{g\infty}\approx(-0.5,0)$ (cells/it.) .

VII.5.2 Stereotype 2

The parameter set (50$\times$50, BHL+L, DAU, OP, 256%, 1, 0%) exhibits the highest transport rate when BHL+L is used and represents another stereotypical attractor state. A snapshot of an attractor is shown in Fig.S4. Such attractors exhibit a transport velocity of $\approx(-0.5,0)$ . The corresponding distribution of the relative power spectral density P$(k_{j},k_{i})$ is shown in Fig.S13. The weighted mean wave vector is directed to the right top and amounts to: $\langle\vec{k}\rangle\approx(0.25,0.25)$.

VII.5.3 Stereotype 3,4

When using BHL+L together with SRAU1 and SRAU2 the model self-organizes towards two further stereotypical stable attractor states. In the accompanying paper one representative snapshot of each type has been illustrated.

Since SRAU1 as well as SRAU2 impose a wave-like activation of the actuators the wave velocity is purely imposed by the activation scheme, which is why we shall skip the discussion considering the characteristics of the exhibited wave field for these two stereotypes.

VII.5.4 Stereotype 5

A small subset of settings generates functioning attractors with toric boundary conditions. A representative is given by the parameter set (50$\times$50, BHL+L, DAU, TO, 256%, 1, 0%). In contrast to the other stereotypes, reached by settings using the prescribed cyclic four-phase motion, these states keep re-organizing and thus, the spatio-temporal correlation function decreases with increasing time lags. The corresponding attractors display many defects (as indicated by the snapshot provided in Fig.S14) and a sequence of states appears rather chaotic, while the propagation of organized wavelets heading towards the right top is clearly recognizable. Accordingly, the relative power spectral density illustrated in Fig.S14 displays a clear peak wave vector: $\vec{k}_{peak}=(0.1,0.25)$. The exhibited terminal transport velocity amounts to $\vec{v}_{g\infty}=(0.3,0)$ and therefore encloses an angle of $\approx\pi/3$ with the peak wave vector.