Evolution of default genetic control mechanisms

William Bains1*, Enrico Borriello2, Dirk Schulze-Makuch3,4,5

1 Department of Earth, Atmospheric and Planetary Sciences, Massachusetts
Institute of Technology. Cambridge. MA 02139, USA

2 ASU-SFI Center for Biosocial Complex Systems, Arizona State University,
Tempe, AZ 85281-2701 , USA

3 Zentrum für Astronomie und Astrophysik, Technische Universität Berlin,
Straße des 17. Juni 135, 10623 Berlin , Germany

4 German Research Centre for Geosciences (GFZ), Section Geomicrobiology,
Potsdam, Germany.

5 Department of Experimental Limnology, Leibniz-Institute of Freshwater
Ecology and Inland Fisheries (IGB), Stechlin, Germany.

* Correspondence to W.Bains<EMAIL_ADDRESS>
## Abstract

We present a model of the evolution of control systems in a genome under
environmental constraints. The model conceptually follows the Jacob and Monod
model of gene control. Genes contain control elements which respond to the
internal state of the cell as well as the environment to control expression of
a coding region. Control and coding regions evolve to maximize a fitness
function between expressed coding sequences and the environment. 118 runs of
the model run to an average of $1.4\times 10^{6}$ ‘generations’ each with a
range of starting parameters probed the conditions under which genomes evolved
a ‘default style’ of control. Unexpectedly, the control logic that evolved was
not significantly correlated to the complexity of the environment. Genetic
logic was strongly correlated with genome complexity and with the fraction of
genes active in the cell at any one time. More complex genomes correlated with
the evolution of genetic controls in which genes were active (‘default on’),
and a low fraction of genes being expressed correlated with a genetic logic in
which genes were biased to being inactive unless positively activated
(‘default off’ logic). We discuss how this might relate to the evolution of
the complex eukaryotic genome, which operates in a ‘default off’ mode.

##  1 Introduction

Obligate multicellularity is uniquely a eukaryotic trait (1-3), and with it
the morphological complexity that comes from combining many distinct cell
types into one organism. Multicellularity requires complex genetic controls
both to provide the control to generate different genetic activity patterns in
different cell types and to provide the ‘programme’ to construct the adult
organism. In addition, the more complex internal architecture and controls in
the eukaryotic cell also require specific controls. In some single-celled
eukaryotes such internal complexity resembles that of equivalently sized
multicellular organisms. Reflecting this, genome sizes in eukaryotes can
exceed those of the largest bacterial or archaeal (“prokaryotic”) genomes by 4
orders of magnitude (Figure 1).

Figure 1: Sizes of completed genome sequences, showing distinct sizes for
prokaryotes (bacteria and archaea), eukaryotes and viruses. $X$ axis: genome
size in megabases. $Y$ axis: fraction of each of the four classes of organism
that have genome of that size. Data from
https://www.ncbi.nlm.nih.gov/genome/browse/#!/overview/, accessed 15th June
2020; based on 27308 bacteria, 1769 Archaea, 5300 Eukarya and 19536 viruses.

There is substantial overlap in coding capacity between the larger prokaryotic
genomes and eukaryotic genomes. The coding capacity of some of the larger
prokaryotic genomes such as those of some cyanobacteria ($\sim 12^{\prime}000$
coding sequences (CDS) (4)), Ktedonobacter racemifer ($\sim 11^{\prime}500$
CDS (5)), Sorangium cellulosum ($\sim 9000$ CDS (6)), Magnetobacterium
bavaricum ($\sim 8500$ CDS (7)) overlaps with coding capacity of multicellular
fungi ($5000-15^{\prime}000$ ( e.g. (8, 9)) and autotrophic protists
($10^{\prime}000-20^{\prime}000$ CDS (9)) and approaches that of Drosophila
melanogaster ($\sim 16^{\prime}000$ CDS (10)). The size difference between
prokaryotic and eukaryotic genomes is primarily due to non-coding DNA that is
related in part to gene control. Thus the E.coli genome has little non-coding
DNA, and $\sim 285$ proteins are involved in gene control (11), $\sim 7\%$ of
the genome. By contrast over 90% of the human genome is non-coding, and
conservative estimates are that 10 times as many non-coding bases as coding
bases are evolutionarily conserved (i.e. are presumed to have selectable
function) (12, 13).

What enabled this increase in genetic complexity? The key difference between
prokaryotic and eukaryotic cells have been suggested to be chemistry,
intracellular structure, energetics and genetics. In general, any small
molecule structure made by a eukaryotic cell will be made by a prokaryote as
well. Many ‘eukaryotic’ cellular structures are actually found in a few
prokaryotes as well. Linear chromosomes are found in bacteria ($14-17$).
Intracellular membrane compartments for secretion and processing (18, 19) and
energy capture (19-24) as well as membrane-bound DNA-containing bodies are
found in Planctomycetes (22, 25). Achromatium oxaliferum contains complex
internal membranes containing calcium carbonate (whose function is obscure)
(26), Entotheonella detoxifies arsenic and barium by sequestering it in
internal vesicles (27), and cyanobacteria have stacked internal photosynthetic
membranes (28). The intracellular membrane system of eukaryotes is integrated
into a dynamic network of vesicle trafficking and control which is rare in
prokaryotes (reviewed in (29)); however some of the core proteins and
structural elements of a cytoskeleton are also found in prokaryotes ($30-35$),
and the giant bacterium Epulopiscium fishelsoni has an internal tubule system
so similar to eukaryotes that it was initially mistaken for a protozoan (36,
37). These examples all suggest that complex structure per se follows from
large size, rather than large size following from internal structure.

It is widely accepted that the modern eukaryotic cell evolved by a series of
endosymbiotic events (38, 39). Recent insights gained from molecular biology
show examples of endosymbiotic bacteria that live inside other bacteria
($40-43$), and bacteria that live inside modern mitochondria (44) (as well as
a wealth of endosymbiotic bacteria in eukaryotic cells) which suggests that
prokaryotic endosymbiotic events, while unusual, are not extremely rare. Lane
and Martin ($45-47$) suggest that the endosymbiotic event, and consequent
development of internal membrane-bound energy-generating organelles, enabled
the ability to generate energy from intracellular membranes acquired through
endosymbiosis is key, as more genes imply more proteins and proteins require
energy to make. We find this theory lacking for three reasons. Firstly, the
majority of genes in the larger eukaryotic genomes do not code protein –
complexity comes from non-coding RNA genes and regulatory elements as
discussed above. Secondly, most of the coding genes in any one cell are not
transcribed; indeed the whole reason to maintain a complex genetic apparatus
is so that different subsets of genes can be expressed at different times.
Genomes containing more coding sequences do not make more proteins at any one
time. Lastly, protein synthesis is only the major use of cellular energy in
autotrophic bacteria grown under conditions of unlimited nutrition. Under more
normal conditions of growth, protein synthesis rarely is observed to consume
more than 20% of the cell’s energy, and of course in non-growing cells (which
is most cells in the biosphere most of the time) protein synthesis is only
needed for maintenance and turnover, a minor part of the overall ‘maintenance
energy’ ($48-50$). (See Supplementary File 1 for a more detailed analysis
protein synthesis’ energy requirements).

We have recently suggested that the default logic of gene control is a
significant factor enabling eukaryogenesis (51). It is observed that it is
easier for a gene to be expressed in a prokaryote than a eukaryote, as
evidenced by the fate of pseudogenes, the construction of expression vectors,
and the fate of differentiated gene expression in cell fusion (reviewed in
(52)) as well as arguments from the mechanisms of gene control (See below).
While there are exceptions, it is broadly true that the eukaryotic genome is
by default ‘off’ and needs metabolic energy to turn ‘on’, whereas in a
prokaryotic genome genes are by default ‘on’ unless turned ‘off’.

This default control mode is reflected in the thermodynamics of gene control.
In eukaryotes, complexes of proteins are required to remodel chromatin around
promoters before genes can be transcribed, involving the ATPase molecular
motors Snf2 and Sth1 (53, 54), and subsequent ATP-dependent binding of
transcription factors to chromatin (55) before RNA polymerase can bind to a
promoter. The remodelling involves (inter alia) ATP-dependent removal of H2A/B
dimers from nucleosomes (56-59) (Archaeal nucleosomes lack H2A/B dimers, and
consist of homologues of H3/4 dimers only (60, 61)). By contrast the molecular
rearrangements that control initiation of bacterial transcription are powered
by the binding energy of the various proteins (11, 62, 63). Archaea have
similar transcription initiation logic to bacteria, despite having RNA
polymerase complexes similar to those in eukaryotes (64, 65). RNA elongation
in bacteria requires roughly 1.5 ATP per base added, again being controlled by
protein binding factors (11). In Eukaryotes ATP-dependent chromatin
remodelling is required for RNA elongation, as well as energy-consuming
histone acetylation and methylation chemistry (66).

We argue that this ‘Default off’ logic is more efficient if the majority of
the genome is silent, as would be the case if the genome encoded many
expression programmes only one of which is active at once. It would also allow
the facile accumulation of silenced duplicate genes to act as the substrate
for genome complexification, which itself is associated with rapid
diversification and adaptation (67). We propose that a ‘default off’ logic
will favour the evolution of complex genomes which code for multiple
expression patterns, a ‘default on’ logic will favour the evolution of
compact, efficient genomes with relatively few distinct phenotypes.

This hypothesis should be testable by simulation and by experiment. As a first
step in this we present a simplified model of gene control and evolution that
can evolve either ‘Default On’ or ‘Default Off’ logic. In this paper we
present the model, and initial results from its execution.

A) 0

B) 0

Figure 2: Summary of model structure. A) overall design philosophy, showing
feedbacks between genotype, its encoded phenotype, and the environment that it
must fit. B) summary of model components. See text for details.

##  2 Methods

###  2.1 Modelling approach

We attempt to model the evolution of control logic of genes under selective
pressure. As a balance between the need for computability on one hand and the
need for biological ‘realism’ on the other, we chose the ‘classical’ operon as
a model on which to build the model structure. A series of sequences upstream
of the coding sequence can bind proteins which allow, promote or catalyse
transcription (positive elements) or which can bind proteins that retard or
prevent transcription (negative elements). A similar process applies to
eukaryotic genes in that positive and negative regulatory elements influence
the transcription of the gene, although in eukaryotes those regulatory
elements may be distant from the gene. Whether those regulatory elements are
active will depend on the proteins in the cell, so that there is feedback
between the phenotype and the transcription of the genotype that it encodes.
The fitness of an organism depends on the ‘fit’ between its phenotype and its
environment, but that environment can change, so the expression of genes must
also be influenced by the environment. The model must also be able to be
queried for some surrogate of ‘default off’ or ‘default on’ genetics
independent of how many genes in an organism are actually transcribed at any
one time (which will depend on the demands of the environment).

The properties of the model are summarised in Figure 2A.

###  2.2 Specifics of the model

To capture the requirements above, the model was constructed as follows. For
simplicity, everything in the model is strings of one type. Thus the phenotype
is a set of strings of the same sort of as the genotype. The strings are made
up of different characters; there can be any number of types of characters (if
the strings were to mimic DNA or RNA, the number of character types would be
4; the model was run with the number of character types ranging from 2 to 16).
There is no equivalent of protein translation in the system. The model
consisted of a number of organisms – in this initial implementation there were
only 5 organisms for computational reasons. The organisms exist in an
environment. Each organism contains a number of genes which together comprise
its genotype; in the runs reported here, organisms contained 25, 50 or 100
genes. Each gene is composed of up to ten positive regulatory elements, up to
ten negative regulatory elements, and a coding sequence. The sum of the coding
regions of genes that are active at any one time comprise the organism’s
phenotype. The organism’s fitness is the match between its phenotype and its
environment as follows. The environment comprises positive elements, negative
elements, and signalling elements. Fitness is the sum of the number of
positive environmental elements that match the current phenotype minus the
number of negative environmental elements that match the current phenotype.
(This is to reflect that sometimes having a function in a cell can be
detrimental to the cell; were this not true, in our model all cells would
express all genes all the time for maximum fitness.)

Gene expression is controlled as follows. A regulatory element is active when
either it matches an environmental signalling element or it matches the
phenotype. This represents the transduction of an environmental signal into
gene activity, and the transduction of internal gene activity into gene
activity. If the sum of the number of active positive regulatory elements
exceed the number of active negative regulatory elements then the gene it
transcribed and its coding sequence is added to the phenotype.

The model is seeded with random strings. At each cycle a new phenotype is
computed, and new fitness computed for each organism , and the most fit
organism randomly replaces one of the other organisms (which can include self-
replacement). The organisms are then mutated by making small, random changes
(character changes, insertions or deletions, with a bias of 6:4 deletion over
insertion) to a fraction (typically between $10^{-5}$ and $5\times 10^{-5}$)
of the strings in the genotype, or completely deleting one of them (typically
with a probability between $10^{-6}$ and $5\times 10^{-6}$).

The model components are summarised in Figure 2B.

##  3 Results

###  3.1 Modelling selection and adaptation

We begin by showing that the model produces results that are consistent with
adaptation, i.e. with changing from an initial random state to a state where
the average fitness of the organisms is greater than it was at the start. We
emphasise that changes made to the components of the model are entirely
random; there is no directionality in the model except a slight bias towards
gene shrinkage noted below. Both the initial genome and the environmental
factors that the genome has to adapt to are randomly generated as well.
Adaptation is therefore the result of selection for better ‘fitness’.

We can measure the ‘degree of perfection’ $P$ of an organism in terms of the
evolved fitness $F$ as a fraction of the possible maximum fitness, as the
maximum fitness is the number of environmental factors $E_{f}$. Some example
fitness curves are shown in Figure 3. Figure 3A shows a typical curve that
reaches a plateau of fitness and then does not achieve any greater fitness in
the run. Figure 3B shows a curve that is similar to $\sim 500^{\prime}000$
generations, but then a new increase in fitness is observed. Figure 3C shows
the decomposition of the fitnesses to each of three environments in a model,
together with the average across all environments. In this run, the organism
is tested against one of three, unrelated environments; the environment that
the organism has to match changes every two generations. Note that fitness to
each environment does not increase in parallel – sometimes selection has
resulted in better fitness for one environment, sometimes for another. Fitness
for an environment can actually decline if overall fitness does not decrease
substantially. Figure 3D shows the separate fitness trajectories of five
organisms as they evolve in a single environment. Again, individual organism
can lose fitness, but the population trend is usually to increasing fitness.
Lastly, Figure 3E shows a model that has not evolved significantly. Most of
the change in fitness in Figure 3E appear to be noise, and fitness wanders
around a low average ($P\sim 0.04$ in this case, as $E_{f}=100$).

A) 0 | B) 0
---|---
|
C) 0 | D) 0
|
E) 0 |
|
Figure 3: Examples of fitness plots for different runs of the model. A)
Average fitness of the population converges smoothly on a maximum. B) Average
fitness shows a jump in adaptation at $500^{\prime}000$ generations. C)
convergence of average fitness across three environments, showing divergent
adaptation to each of the environments. D) Fitness of each of the five
organisms making up the population plotted separately in a run that converges
on a solution. E) Plot of fitness in a run that fails to converge on an
optimum fitness.

###  3.2 Failure to adapt

Models did not converge onto a fit state $\sim 1/5$ of the time (depending on
what ‘fit’ means, and the selection of parameters). This was found to be a
function of the degree to which the genome complexity can match the
environmental complexity (Figure 4). Highly complex environments are not
efficiently matched by low complexity genomes, as would be expected. This to a
degree is a consequence of limited run-time: some model runs had low fitness
for a time and then experienced a ‘jump’ in fitness as a low-probability
solution was ‘discovered’ (e.g. Figure 3B).

We define whether a population is converging on a solution with a Curve
Parameter $C_{p}$ as follows: we define a time $P$ as the time at which the
population reaches a plateau of adaptation, i.e. does not appear by inspection
to be able to increase its adaptation. $C_{p}$ distinguishes between
populations that smoothly approach such a fitness plateau, such as shown in
Figure 3A, and populations whose fitness fluctuates, such as in Figure 3E.
Thus, if the fitness at times $0.25P$, $0.5P$, $0.75P$ and $P$ are $A$, $B$,
$C$, $D$ respectively, and $S(x)$ is the sign of $x$, (such that
$x>0\Rightarrow s(x)=1$; $X<0\Rightarrow s(x)=-1$; $x=0\Rightarrow s(x)=0$)
then

$C_{p}=s(B-A)+s(C-A)+s(D-A)+s(C-B)+s(D-B)+s(D-C).$

If $A<B<C<D$ (i.e. fitness is increasing throughout the run), then $C_{p}=6$.

For this analysis, runs of the model were only used if the curve parameter
$C_{p}$ is greater than zero. 101 out of 118 runs of the model met this
criterion. Omitting curves for which $C_{p}\leq 0$ resulted in substantially
less scatter in the adaptation curves averaged across all models (Figure 5).

Would all models converge on an optimal solution eventually? We hypothesise
they would, but it might take years to achieve this using the relatively
inefficient coding, which for practical reasons was run for only an average of
$1.4\times 10^{6}$ generations. (For comparison, the long-term evolution
experiments performed by the Lenski lab. have been running for more than
$20^{\prime}000$ generations, and show a range of adaptations in gene control
structure without changing the underlying mechanisms or logic of the gene
control architecture ($68-70$)). Models were therefore stopped when they
seemed to have reached a steady state of control structure (as defined below)
and fitness within this time window. Future work will seek a more objective
measure of termination state through spectral analysis of the fitness
functions in an exhaustive scan of our parameter space, Disentangling the
typical timescale of purely stochastic fluctuations from the timescale of
selection-induced changes will provide a useful estimate of the average,
expected time until convergence, as well as estimates of the range on that
time.

Figure 4: Selection is inefficient for runs with a combination of high
environmental complexity and low genome complexity. $X$ axis: genome
complexity (number of genes in each organism times the number of types of
bases of which those genes are made up, times the average length of the genes
at the end of the selection process. $Y$ axis: environmental complexity;
number of environmental factors that must match the genotype, multiplied by
the number of different environments to be fitted, the number of base types in
each sequence, and the length of the environmental strings to be matched. The
‘perfection index’ –fitness at selection plateau divided by maximum possible
fitness– is both proportional to the circle areas and, for enhanced
readability, to the colour scale (vertical bar).  Figure 5: Average fitness of
the model runs. Curves are normalized to maximum fitness = 1. The $X$ axis
shows the fraction of time until the fitness reaches a stable plateau. The
average fitness across all 118 model runs was calculated for time
$\in[0,0.25]$, time $\in[0.25,0.5]$ , time $\in[0.5,0.75]$, time $\in[0.75,1]$
and $\textrm{time}>1$ (i.e. on the plateau of fitness). By definition, fitness
= 1 when $\textrm{time}>1$. Error bars are standard deviations. Blue curve:
All 118 model runs. Orange curve: 101 model runs for which Curve Parameter
$C_{p}>0$.

###  3.3 Coding region evolution

There is a slight bias the mutation mechanism towards gene shrinkage, included
because a) this is seen in real mutation rates and b) it protects the model
against indefinite expansion of genes through ‘drift’. Despite this, the
average length of coding regions tends to increase with model progression
(Figure 6). This is explicable as follows. Gene activation depends on matching
part of an expressed coding sequence to a regulatory sequence. Thus larger
genes mean a greater chance of productive interaction with a regulatory
element. The only selective pressure against long genes is the chance that
they interact with one of the ‘negative’ environmental elements when they are
expressed. As for all runs there are more regulatory elements (20 per gene)
than environmental factors (a maximum of 200), this provides a selection
pressure towards longer genes.

Figure 6: Length of the coding sequences averaged across all genes in a run
($Y$ axis) as a function of the starting length of those sequences ($X$ axis).
Circle sizes and color scale show the number of genes in the run, and show no
distinct pattern.

###  3.4 Default genetic control measures

The computational effort to exhaustively analyse the entire control network of
up to 100 genes each interacting with up to 20 control elements in each gene
in each of 5 organisms in 100 runs or up to $50^{\prime}000$ time steps each
is unrealistic, and so we summarise the overall style of control as follows.

Average control element length. A short regulatory element is more likely to
match a sequence in the phenotype than a long regulatory element, because a
short string is more likely to match a random target by chance than is a long
string. (Consider the chance that the strings “A” and “ALPHABET” will match
the text in this paper) Thus if negative regulatory elements are on average
shorter than positive regulatory elements in a gene, it is likely that the
gene is not active. Thus the ratio

$R_{a}=\textrm{ ratio (average regulatory element length)}$ $\
=\frac{\sum\textrm{neg. elem. len.}/\sum\textrm{neg. elem.
count}}{\sum\textrm{pos. elem. len.}/\sum\textrm{pos. elem. count}}$

is a measure of the bias towards inactive genes, i.e. of ‘default off’
regulatory style. (Zero-length elements, i.e. ones which have been deleted,
are not counted in the average).

Average minimum control element length. The problem with average regulatory
element length as a measure of control logic is that one short regulatory
element (likely to be active) can dominate the control of a gene over a number
of long regulatory elements (which are in effect ‘junk DNA’, never being
active). So the shortest regulatory element is the one most likely to be
‘active’ in a gene. If the shortest positive regulatory element is shorter
than the shortest negative regulatory element, then there is a greater chance
that the gene will be active; if the shortest negative regulatory element is
shorter than the shortest positive element, then the gene is more likely to be
inactive. We therefore adopted a measure looking for the shortest regulatory
element in a gene. For each gene, the shortest non-zero control sequence is
recorded for positive and negative control elements. The average of the length
of the shortest regulatory element for all genes in the genotype is reported.
Thus the ratio

$R_{m}=\textrm{ ratio (minimum regulatory element length)}$ $\
=\frac{(\sum\textrm{min. neg. elem.
len./gene})/\sum\textrm{genes}}{(\sum\textrm{min. pos. elem.
len./gene})/\sum\textrm{genes}}$

is a measure of the bias towards inactive genes, i.e. of ‘default off’
regulatory style. (Again, zero-length elements, i.e. ones which have been
deleted, are not counted in the average)

Correlations of these two measures to both the inputs and the outputs of model
runs are provided in Table 1. We emphasise that this is an initial modelling
study, and much more extensive modelling with more efficiently coded models
and better hardware will be needed to confirm, expand and dissect these
findings. However three patterns are clear from the results here.

Model parameters | Measures of genetic control logic
---|---
|  | min $-$ve / min $+$ve
---
($<1=$ default off)
| avg -ve / avg $+$ve
---
($<1=$ default off)
Starting parameters
Number of environments |  $-$0.02 | 0.032
Number of environmental factors | 0.262 * | 0.144
| Environmental complexity
---
(a)
0.141 | 0.151
Length of initial gene | 0.445 ** |  $-$0.044
Number of genes |  $-$0.144 |  $-$0.144
Genome complexity at start (b) | 0.297 ** | 0.096
Parameters at fitness plateau
Adaptation at plateau | 0.326 ** | 0.167
Fraction of perfection | 0.117 |  $-$0.015
Average gene length | 0.683 **** | 0.198
Genome complexity at end (b) | 0.373 ** | 0.180
regulatory complexity at end (c) | 0.13 | 0.184
Number of genes expressed | 0.262 * | 0.102
Fraction of genes expressed | 0.713 **** | 0.454 **
Table 1: Correlations with control logic style.

Correlations between two measures of genetic ‘style’ and some inputs and
outputs from models.

  * (a)

Environmental complexity is the [number of characters] $\times$ [number of
environments] $\times$ [number of factors per environment];

  * (b)

Genome complexity = [number of characters] $\times$ [number of genes] $\times$
[average length of coding regions];

  * (c)

Regulatory complexity = [number of characters] $\times$ [number of genes]
$\times$ [average length of regulatory elements].

“Min $–$ve” = average length of the shortest negative regulatory element in
each gene, averaged over all genes. “min $+$ve“= average length of the
shortest positive regulatory element in each gene, averaged over all genes.
“Avg $–$ve” = average length of all non-zero-length negative regulatory
elements in the genome.”Avg $+$ve”= average length of all non-zero-length
positive regulatory elements in the genome. Significance of the correlations
(i.e. chance that the observed correlation is seen in 101 model runs if the
measures of control logic are not correlated with the input parameters)
$*=p<0.05,**=p<0.01,***=p<0.001,***\,*=p<0.0001$.

Firstly, whether default on or default off genetics evolves is largely
independent of how complex the environment is. This was a surprising result,
as we expected more complex environments to drive selection for more complex
genetic controls, with consequences (positive or negative) for a ‘default off’
control style. From this initial analysis, such effects do not appear to
dominate. This may be an artefact of the small sizes of the populations
(leading to noise which obscures patterns), the small number of combinations
of parameters explored (61 out of 720 possible combinations of the parameters
used in these runs), or the small genomes (maximum 100 genes).

Secondly, genome complexity both at the start and at the end of the models is
correlated with ‘default on’ control style. Again, this was a surprise. Our
hypothesis is that ‘default off’ genetics allows complex genome evolution.
However, our initial hypothesis, that ‘default off’ genetics allows ready gene
duplication, is not captured in this model, where the number of genes is
fixed.

Lastly, ‘default on’ genetics is most strongly correlated with the number of
expressed genes. We further dissect this in Figure 7. There is a striking
correlation between two measures of ‘default control logic’ and the number of
expressed genes. If relatively few genes are expressed in a genome, then
‘default off’ is preferred. If many genes are expressed, ‘default on’ is
preferred. Note that the correlation with the number of expressed genes is
much weaker (Table 1)– it is the fraction of the genome which is expressed
that correlates more strongly with genetic logic than any other parameter.

Figure 7: Relationship between the ratio of minimum negative elements to
minimum positive elements ($X\textrm{ axis:}<1=$‘default off’) to the fraction
of genes in a genome expressed at fitness plateau ($Y$ axis). Both circle area
and color scale are proportional to genome size (25, 50 or 100 genes). “Min
$–$ve” = average length of the shortest negative regulatory element in each
gene, averaged over all genes. “min $+$ve“= average length of the shortest
positive regulatory element in each gene, averaged over all genes.

###  3.5 Reproducibility

No stochastic model will give the same results in different runs, so it is
important to show that the variability in output is not so extreme as to
render results uninterpretable. The purpose of this modelling was to test the
model concept and provide an initial exploration of parameter space: as a
result, only a few sets of runs of the model were replicate runs with the same
parameters. We chose three runs that gave different control logic outputs in
an initial run and re-ran the same parameters with different starting genomes
and environments. The results are summarised in Figure 8. This shows that,
while results are variable, the genetic control outputs, and specifically the
$R_{m}$ parameter, are consistent within replicates: Replicates of a model run
that gave $R_{m}=1$ consistently gave $R_{m}\sim 1$ (“Replicate 1”),
Replicates of a run which gave $R_{m}>1$ (“Replicate 2”) consistently gave
$R_{m}>1$, and Replicates of a run that gave $R_{m}<1$ consistently gave
$R_{m}<1$ (“Replicate 3”)

A) 0

B) 0

Figure 8: Reproducibility across runs. A) Example outcomes from all runs with
diverse starting conditions, and B) From replicate sets of runs started from
the same set of parameters. $X$ axis: ‘Perfection index’ (fitness at the
fitness plateau as a fraction of the maximum possible fitness with those
parameters). $Y$ axis: ratio of the length of the average minimum negative
regulatory element length to the average minimum positive regulatory element
length. Both circle area and color scale are proportional to the number of
steps taken to reach a fitness plateau.

##  4 Discussion and Conclusions

We have presented a model of the evolution of genome control logic, and an
initial analysis of its performance on a small number of test cases. The model
performs in a comprehensible way, and evolves fitter organisms. Preliminary
statistics suggest that the model performance is stable, i.e. a given set of
starting conditions will give a set of outputs more closely related to each
other than random, despite the model being a stochastic one.

We emphasise that this is a preliminary exploration of this model only, and
much more needs to be done. However with that caveat, the results show three
things of potential interest to the hypothesis that stimulated its creation

  * i

The genetic logic a population of organisms evolves is not related to the
complexity of the environment it finds itself in. This was unexpected

  * ii)

The evolved genetic logic is strongly related to the starting and the final,
evolved genome complexity. More complex genomes have ‘default on’ logic. This
is not predicted by the model, but as the model’s predictions on evolution of
genetic logic refers primarily to the acquisition of new genes in the genome,
an aspect of evolution not captured here, this does not test the hypothesis.

  * iii)

The strongest correlations with genetic logic are with the fraction of the
genome that is expressed.

Point (iii) above fits with (although is a weak test of) our original
hypothesis. It also fills in a significant gap in the hypothesis about why a
‘default off’ logic should be selected. Clearly, an organism cannot evolve
‘default off’ in anticipation of acquiring new genes. However if a specific
combination of environmental and genetic features encouraged the development
of ‘default off’ genetics, then such an organism would be pre-adapted for
genome complexification by gene duplication and divergence. As noted in the
introduction, the majority of genes in eukaryotic genomes are not expressed at
any one time. Most of them are ‘off’. Our model appears to be evolving a
similar expression pattern in some cases, and in those cases the ‘default off’
logic is selected.

If our results represent the more complex world of real genetics, then we
might speculate that organisms living in an environment that occasionally
called on a diverse set of genes but most of the time did not require them
would feel short-term selective pressure to evolve a ‘default off’ logic. Such
an environment could be one in which a heterotroph lived in a community made
up of a changing composition of autotrophs, each of which provided a small
number of substrates to the heterotroph. If such a scenario were valid, then
we would expect more comprehensive modelling to reveal influences of
environmental change on both expressed gene numbers and default logic. Such
work is being actively pursued.

## References

  1. 1.

Knoll AH. The Multiple Origins of Complex Multicellularity. Annual Review of
Earth and Planetary Sciences. 2011;39(1):217-39.

  2. 2.

Parfrey LW, Lahr DJG. Multicellularity arose several times in the evolution of
eukaryotes (Response to DOI 10.1002/bies.201100187). BioEssays.
2013;35(4):339-47.

  3. 3.

Rokas A. The Origins of Multicellularity and the Early History of the Genetic
Toolkit For Animal Development. Annual Review of Genetics. 2008;42(1):235-51.

  4. 4.

Dagan T, Roettger M, Stucken K, Landan G, Koch R, Major P, et al. Genomes of
Stigonematalean Cyanobacteria (Subsection V) and the Evolution of Oxygenic
Photosynthesis from Prokaryotes to Plastids. Genome biology and evolution.
2013;5(1):31-44.

  5. 5.

Chang Y-j, Land M, Hauser L, Chertkov O, Glavina Del Rio T, Nolan M, et al.
Non-contiguous finished genome sequence and contextual data of the filamentous
soil bacterium Ktedonobacter racemifer type strain (SOSP1-21T). Standards in
Genomic Sciences. 2011;5(1):97-111.

  6. 6.

Schneiker S, Perlova O, Kaiser O, Gerth K, Alici A, Altmeyer MO, et al.
Complete genome sequence of the myxobacterium Sorangium cellulosum. Nat
Biotech. 2007;25(11):1281-9.

  7. 7.

Kolinko S, Richter M, Glöckner F-O, Brachmann A, Schüler D. Single-cell
genomics of uncultivated deep-branching magnetotactic bacteria reveals a
conserved set of magnetosome genes. Environmental Microbiology.
2016;18(1):21-37.

  8. 8.

Xu J, Saunders CW, Hu P, Grant RA, Boekhout T, Kuramae EE, et al. Dandruff-
associated Malassezia genomes reveal convergent and divergent virulence traits
shared with plant and human fungal pathogens. Proceedings of the National
Academy of Sciences. 2007;104(47):18730-5.

  9. 9.

Anantharaman V, Iyer LM, Aravind L. Comparative Genomics of Protists: New
Insights into the Evolution of Eukaryotic Signal Transduction and Gene
Regulation. Annual review of microbiology. 2007;61(1): 453-75.

  10. 10.

Adams MD, Celniker SE, Holt RA, Evans CA, Gocayne JD, Amanatides PG, et al.
The Genome Sequence of Drosophila melanogaster. Science.
2000;287(5461):2185-95.

  11. 11.

Borukhov S, Lee J, Laptenko O. Bacterial transcription elongation factors: new
insights into molecular mechanism of action. Molecular Microbiology.
2005;55(5):1315-24.

  12. 12.

Lander ES, Linton LM, Birren B, Nusbaum C, Zody MC, Baldwin J, et al. Initial
sequencing and analysis of the human genome. Nature. 2001;409:860.

  13. 13.

Ponting CP, Hardison RC. What fraction of the human genome is functional?
Genome Research. 2011;21(11):1769-76.

  14. 14.

Ferdows MS, Barbour AG. Megabase-sized linear DNA in the bacterium Borrelia
burgdorferi, the Lyme disease agent. Proceedings of the National Academy of
Sciences. 1989;86(15):5969-73.

  15. 15.

Hinnebusch J, Tilly K. Linear plasmids and chromosomes in bacteria. Molecular
Microbiology. 1993;10(5):917-22.

  16. 16.

Sherman L, Min H, Toepel J, Pakrasi H. Better Living Through Cyanothece –
Unicellular Diazotrophic Cyanobacteria with Highly Versatile Metabolic
Systems. In: Hallenbeck PC, editor. Recent Advances in Phototrophic
Prokaryotes. Advances in Experimental Medicine and Biology. 675: Springer New
York; 2010. p. 275-90.

  17. 17.

Kube M, Schneider B, Kuhl H, Dandekar T, Heitmann K, Migdoll A, et al. The
linear chromosome of the plant-pathogenic mycoplasma ’Candidatus Phytoplasma
mali’. BMC Genomics. 2008;9(1):306.

  18. 18.

Kulp A, Kuehn MJ. Biological Functions and Biogenesis of Secreted Bacterial
Outer Membrane Vesicles. Annual review of microbiology. 2010;64:163-84.

  19. 19.

Stolz JF. Bacterial Intracellular Membranes. eLS: John Wiley & Sons, Ltd;
2001.

  20. 20.

Hanson RS, Hanson TE. Methanotrophic bacteria. Microbiological Reviews.
1996;60(2):439-71.

  21. 21.

Prust C, Hoffmeister M, Liesegang H, Wiezer A, Fricke WF, Ehrenreich A, et al.
Complete genome sequence of the acetic acid bacterium Gluconobacter oxydans.
Nat Biotech. 2005;23(2):195-200.

  22. 22.

Fuerst JA. Intracellular compartmentalization in Planctomycetes. Ann Rev
Microbiology. 2005;59:299 - 328.

  23. 23.

van Niftrik LA, Fuerst JA, Damsté JSS, Kuenen JG, Jetten MSM, Strous M. The
anammoxosome: an intracytoplasmic compartment in anammox bacteria. FEMS
Microbiology Letters. 2004;233(1):7-13.

  24. 24.

van Niftrik L, Geerts WJC, van Donselaar EG, Humbel BM, Yakushevska A,
Verkleij AJ, et al. Combined structural and chemical analysis of the
anammoxosome: A membrane-bounded intracytoplasmic compartment in anammox
bacteria. Journal of Structural Biology. 2008;161(3):401-10.

  25. 25.

Fuerst JA, Webb RI, Garson MJ, Hardy L, Reiswig HM. Membrane-bounded nucleoids
in microbial symbionts of marine sponges. FEMS Microbiology Letters.
1998;166(1):29-34.

  26. 26.

Schorn S, Salman-Carvalho V, Littmann S, Ionescu D, Grossart H-P, Cypionka H.
Cell Architecture of the Giant Sulfur Bacterium Achromatium oxaliferum: Extra-
cytoplasmic Localization of Calcium Carbonate Bodies. FEMS Microbiology
Ecology. 2019;96(2).

  27. 27.

Keren R, Mayzel B, Lavy A, Polishchuk I, Levy D, Fakra SC, et al. Sponge-
associated bacteria mineralize arsenic and barium on intracellular vesicles.
Nature Communications. 2017;8(1):14393.

  28. 28.

Kunkel DD. Thylakoid centers: Structures associated with the cyano- bacterial
photosynthetic membrane system. Archives of Microbiology. 1982;133(2):97-9.

  29. 29.

Jékely G. Origin of eukaryotic endomembranes. In: Jékely G, editor. Eukaryotic
Membranes and Cytoskeleton: Origins and Evolution. Berlin: Springer; 2007. p.
38 - 51.

  30. 30.

Roeben A, Kofler C, Nagy I, Nickell S, Ulrich Hartl F, Bracher A. Crystal
Structure of an Archaeal Actin Homolog. Journal of Molecular Biology.
2006;358(1):145-56.

  31. 31.

Vollmer W. The prokaryotic cytoskeleton: a putative target for inhibitors and
antibiotics? Applied Microbiology and Biotechnology. 2006;73(1):37-47.

  32. 32.

van den Ent F, Amos LA, Lowe J. Prokaryotic origin of the actin cytoskeleton.
Nature. 2001;413(6851):39-44.

  33. 33.

Erickson HP. FtsZ, a tubulin homologue in prokaryote cell division. Trends in
Cell Biology. 1997;7(9):362-7.

  34. 34.

Ausmees N, Kuhn JR, Jacobs-Wagner C. The Bacterial Cytoskeleton: An
Intermediate Filament-Like Function in Cell Shape. Cell. 2003;115(6):705-13.

  35. 35.

Williams TA, Foster PG, Cox CJ, Embley TM. An archeal origin of eukaryotes
supports only two primary domains of life. Nature. 2014;504:231 - 6.

  36. 36.

Montgomery WL, Pollak PE. Epulopiscium fishelsoni N. G., N. Sp., a Protist of
Uncertain Taxonomic Affinities from the Gut of an Herbivorous Reef Fish.
Eukaryotic microbiology. 1988;35(4):565 - 9.

  37. 37.

Angert ER, Clements KD, Pace NR. The largest bacterium. Nature.
1993;362(6417):239-41.

  38. 38.

de Duve C. The origin of eukaryotes: a reappraisal. Nature Reviews Genetics.
2007;8:395 - 403.

  39. 39.

Pittis AA, Gabaldón T. Late acquisition of mitochondria by a host with
chimaeric prokaryotic ancestry. Nature. 2016;531(7592):101 - 4.

  40. 40.

Thao ML, Gullan PJ, Baumann P. Secondary ($\gamma$-Proteobacteria)
Endosymbionts Infect the Primary ($\beta$-Proteobacteria) Endosymbionts of
Mealybugs Multiple Times and Coevolve with Their Hosts. Applied and
Environmental Microbiology. 2002;68(7):3190-7.

  41. 41.

Wujek DE. Intracellular bacteria in the blue-green alga Pleurocapsa minor.
Transactions of the American Microscopic Society. 1979;98(1): 143 - 5.

  42. 42.

Larkin JM, Henk MC. Filamentous sulfide-oxidizing bacteria at hydrocarbon
seeps of the Gulf of Mexico. Microscopy research and technique.
1996;33(1):23-31.

  43. 43.

von Dohlen CD, Kohler S, Alsop ST, McManus WR. Mealybug [beta]-proteobacterial
endosymbionts contain [gamma]-proteobacterial symbionts. Nature.
2001;412(6845):433-6.

  44. 44.

Sassera D, Beninati T, Bandi C, Bouman EAP, Sacchi L, Fabbi M, et al.
‘Candidatus Midichloria mitochondrii’, an endosymbiont of the tick Ixodes
ricinus with a unique intramitochondrial lifestyle. International Journal of
Systematic and Evolutionary Microbiology. 2006;56(11): 2535-40.

  45. 45.

Lane N. Energetics and genetics across the prokaryote-eukaryote divide. Biol
Direct. 2011;6:35.

  46. 46.

Lane N, Martin W. The energetics of genome complexity. Nature.
2010;467(7318):929-34.

  47. 47.

Lane N, Martin WF. Eukaryotes really are special, and mitochondria are why.
Proceedings of the National Academy of Sciences. 2015; 112(35):E4823.

  48. 48.

Tijhuis L, Van Loosdrecht MCM, Heijnen JJ. A thermodynamically based
correlation for maintenance gibbs energy requirements in aerobic and anaerobic
chemotrophic growth. Biotechnology and Bioengineering. 1993;42(4):509-19.

  49. 49.

Hoehler T. Biological energy requirements as quantitative boundary conditions
for life in the subsurface. Geobiology. 2004;2(4):205-15.

  50. 50.

Hoehler TM, Amend J, Shock EL. A “Follow the Energy” Approach for
Astrobiology. Astrobiology. 2007;7(6):819-23.

  51. 51.

Bains W, Schulze-Makuch D. Mechanisms of Evolutionary Innovation Point to
Genetic Control Logic as the Key Difference Between Prokaryotes and
Eukaryotes. Journal of Molecular Evolution. 2015:1-20 (DOI
10.1007/s00239-015-9688-6).

  52. 52.

Bains W, Schulze-Makuch D. The Cosmic Zoo: The (Near) Inevitability of the
Evolution of Complex, Macroscopic Life. Life. 2016;6(3):25.

  53. 53.

Wagner FR, Dienemann C, Wang H, Stützer A, Tegunov D, Urlaub H, et al.
Structure of SWI/SNF chromatin remodeller RSC bound to a nucleosome. Nature.
2020;579(7799):448-51.

  54. 54.

Ye Y, Wu H, Chen K, Clapier CR, Verma N, Zhang W, et al. Structure of the RSC
complex bound to the nucleosome. Science. 2019;366(6467): 838-43.

  55. 55.

Cairns BR. Chromatin remodeling machines: similar motors, ulterior motives.
Trends in Biochemical Sciences. 1998;23(1):20-5.

  56. 56.

Kireeva ML, Walter W, Tchernajenko V, Bondarenko V, Kashlev M, Studitsky VM.
Nucleosome Remodeling Induced by RNA Polymerase II: Loss of the H2A/H2B Dimer
during Transcription. Molecular Cell. 2002;9(3):541-52.

  57. 57.

Mizuguchi G, Shen X, Landry J, Wu W-H, Sen S, Wu C. ATP-Driven Exchange of
Histone H2AZ Variant Catalyzed by SWR1 Chromatin Remodeling Complex. Science.
2004;303(5656):343-8.

  58. 58.

Shen X, Mizuguchi G, Hamiche A, Wu C. A chromatin remodelling complex involved
in transcription and DNA processing. Nature. 2000;406(6795):541-4.

  59. 59.

Olave IA, Peck-Peterson SI, Crabtree GR. Nuclear actin and actin-related
proteins in chromatin remodelling. Annual Review of Biochemistry. 2002;71:755
- 81.

  60. 60.

Grayling RA, Sandman K, Reeve JN. Histones and chromatin structure in
hyperthermophilic Archaea. FEMS Microbiology Reviews. 1996;18(2-3):203-13.

  61. 61.

Sandman K, Reeve JN. Archaeal histones and the origin of the histone fold.
Current Opinion in Microbiology. 2006;9(5):520-5.

  62. 62.

Saecker RM, Record MT, deHaseth PL. Mechanism of Bacterial Transcription
Initiation: RNA Polymerase - Promoter Binding, Isomerization to Initiation-
Competent Open Complexes, and Initiation of RNA Synthesis. Journal of
Molecular Biology. 2011;412(5):754-71.

  63. 63.

Browning DF, Busby SJW. The regulation of bacterial transcription initiation.
Nature Reviews Microbiology. 2004;2(1):57-65.

  64. 64.

Reeve JN. Archaeal chromatin and transcription. Molecular Microbiology.
2003;48(3):587-98.

  65. 65.

Geiduschek EP, Ouhammouch M. Archaeal transcription and its regulators.
Molecular Microbiology. 2005;56(6):1397-407.

  66. 66.

Tsukiyama T, Wu C. Chromatin remodeling and transcription. Current Opinion in
Genetics & Development. 1997;7(2):182-91.

  67. 67.

McGee MD, Borstein SR, Meier JI, Marques DA, Mwaiko S, Taabu A, et al. The
ecological and genomic basis of explosive adaptive radiation. Nature.
2020;586(7827):75-9.

  68. 68.

Crozat E, Philippe N, Lenski RE, Geiselmann J, Schneider D. Long-term
experimental evolution in Escherichia coli. XII. DNA topology as a key target
of selection. Genetics. 2005;169(2):523-32.

  69. 69.

Philippe N, Crozat E, Lenski RE, Schneider D. Evolution of global regulatory
networks during a long-term experiment with Escherichia coli. Bioessays.
2007;29(9):846-60.

  70. 70.

Lamrabet O, Plumbridge J, Martin M, Lenski RE, Schneider D, Hindré T.
Plasticity of promoter-core sequences allows bacteria to compensate for the
loss of a key global regulatory gene. Molecular biology and evolution.
2019;36(6):1121-33.


# Segregation Dynamics with Reinforcement Learning and Agent Based Modeling

Egemen Sert  New England Complex Systems Institute, Cambridge, MA  Department
of Electrical and Electronics Engineering, Middle East Technical University,
Ankara, Turkey  Yaneer Bar-Yam  New England Complex Systems Institute,
Cambridge, MA  Alfredo J. Morales  New England Complex Systems Institute,
Cambridge, MA  MIT Media Lab, Cambridge, MA.  Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

Societies are complex. Properties of social systems can be explained by the
interplay and weaving of individual actions. Incentives are key to understand
people’s choices and decisions. For instance, individual preferences of where
to live may lead to the emergence of social segregation. In this paper, we
combine Reinforcement Learning (RL) with Agent Based Models (ABM) in order to
address the self-organizing dynamics of social segregation and explore the
space of possibilities that emerge from considering different types of
incentives. Our model promotes the creation of interdependencies and
interactions among multiple agents of two different kinds that want to
segregate from each other. For this purpose, agents use Deep Q-Networks to
make decisions based on the rules of the Schelling Segregation model and the
Predator-Prey model. Despite the segregation incentive, our experiments show
that spatial integration can be achieved by establishing interdependencies
among agents of different kinds. They also reveal that segregated areas are
more probable to host older people than diverse areas, which attract younger
ones. Through this work, we show that the combination of RL and ABMs can
create an artificial environment for policy makers to observe potential and
existing behaviors associated to incentives.

##  1 Introduction

The recent availability of large datasets collected from various resources,
such as digital transactions, location data and government census, is
transforming the ways we study and understand social systems [1]. Researchers
and policy makers are able to observe and model social interactions and
dynamics in great detail, including the structure of friendship networks [2],
the behavior of cities [3], politically polarized societies [4], or the spread
of information on social media [5]. These studies show the behaviors present
in the data but do not explore the space of possibilities that human dynamics
may evolve to. Robust policies should consider mechanisms to respond to every
type of events [6], including those that are very rare [7]. Therefore it is
crucial to develop simulation environments such that potentially unobserved
social dynamics can be assessed empirically.

Agent Based Modeling (ABM) is a generative approach to study social phenomena
based on the interaction of individuals [8]. These models show how different
types of individual behavior give rise to emergent macroscopic regularities
[9], such as unequal wealth distributions [10], new political actors [11],
multipolarity in interstate systems [12] and cultural differentiation [13].
Moreover, ABM allows testing core sociological theories against simulations
[10] with emphasis on heterogeneous, autonomous actors with bounded, spatial
information [14]. However, the rules of agent interactions are generally fixed
which limits the exploration of the space of possible behaviors.

Reinforcement Learning (RL) is a simulation method where agents become
intelligent and create new, optimal behaviors based on the state of their
environment and a previously defined structure of incentives. This method is
referred as Multi-Agent Reinforcement Learning (MARL) if multiple agents are
employed. Recently, the combination of RL with Deep Learning architectures
achieve human level performance in complex tasks, including video gaming [15],
motion in harsh environments [16], and effective communication networks
without assumptions [17]. Moreover, it has been recently applied to study
societal dilemma and game theory problems [18] such as the emergence of
cooperation [19, 20], the Prisoner’s Dilemma [21] and payoff matrices in
equilibrium [22]. Although Deep RL algorithms applied to multiple agents
(MARL) can shed light on social phenomena, to the best of our knowledge, the
applications of these methods has been confined to classical game-theoretic
problems [23] and drawing connections to real-world examples remains
unexplored.

In this paper we extend the standard ABM of social segregation using MARL in
order to explore the space of possible behaviors as we modify the structure of
incentives and promote the interaction among agents of different kinds. The
idea is to observe the behavior of agents that want to segregate from each
other when interactions are promoted. We achieve the segregation dynamics by
considering the rules from the Schelling model [9]. The creation of
interdependencies among agents of different kinds is inspired by the dynamics
of the Predator-Prey model [8] where agents hunt each other. Our experiments
show that spatial segregation diminishes as more interdependencies among
agents of different kinds are added. Moreover, our results shed light on
previously unknown behaviors regarding segregation and the age of individuals
which we confirmed using Census data. These methods can be extended to study
other type of social phenomena and inform policy makers on possible actions.

The organization of the paper is as follows: In Section 2 we explain the
experimental setup. Section 3 illustrates the experiment outcomes. In Section
4 we conclude and discuss our results. Future improvements are presented in
Section S1 in the Supplement.

##  2 Methods

Figure 1: Schematic of the model simulation and network architecture. Top
panel: Grid world of experiments. The grid size is 50x50 locations. Red and
blue squares denote two types of agents respectively. White represents empty
regions. Each type has its own Deep Q-Network. Every agent has a field of view
of 11x11 locations. Green border denotes the field of view of the agent
illustrated in green. Agents can move across empty spaces. Bottom panel:
Network structure of $\phi_{A}$ and $\phi_{B}$. Each network receives an input
of 11x11 locations, runs it through five convolution steps and concatenates
the resulting activations with the agent’s remaining age normalized by the
maximum initial age. The feature vector is mapped over the action space using
a fully connected layer. The action with the maximum Q-value is taken for the
agent. Table 1: Training parameters of the Deep Q-Networks used during the
experiments.
Parameter | Value |
---|---|---
Number of Episodes | 1 |
Batch Size | 256 |
Number of Iterations | 5000 |
Number of Training Steps | 60.000 |
Experience Memory Length | 1.000.000 |
Discount Factor ($\gamma$) | 0.98 |
Learning Rate | 0.001 |
Momentum | 0.999 |
Double Network Copy Parameter ($\tau$) | 0.05 |
Initial Exploration Rate | 0.999 |
Final Exploration Rate | 0 |
Exploration Decay (per agent action) | 100.000 |

We design a game in which agents are promoted to both self-segregate and
interact with others. By varying the reward of interactions we are able to
explore different incentives that affect the self-organizing process of
segregation. Our experiments are based on two types of agents: A and B. Agents
try to survive in a 50x50 grid where they can move around and interact with
other agents. They observe an 11x11 patch of the grid centered around their
current position and can live for a total of 100 iterations in isolation.
Figure 1 shows an schematic view of the grid world and the agents. Distinct
colors indicate the agents’ types and the green square represents the
observation window of the agent illustrated in green.

Figure 2: Agents collective behavior for multiple values of interdependence
reward (rows) at multiple times (columns). Rows represent outcomes associated
to different values of the interdependence reward (IR). Columns show the state
of the system at different points of the simulation. Experiments are
initialized with equal initial conditions and random seed. The heat maps are
obtained by averaging over the last 1000 iterations. In Panel (a) red regions
denote biased occupation of type A agents where areas fully occupied with type
A agents are indicated by type concentration of +1. Blue regions denote biased
occupation of type B agents where full occupation of blues are indicated by
type concentration of -1. White areas indicate uniform mixing across types,
indicated by type concentration of 0. In Panel (b) color indicates the age of
agents irrespective of their type. As color shades from blue to red agent age
increases. In Section 3, we introduce a set of videos that represent the
experiments used in creating the heat maps. Figure 3: Segregation dynamics for
multiple values of interdependence reward (IR). Colors correspond to the
results for multiple values of interdependence reward, ranging from yellow
(low) to black (high). The curves are obtained by averaging 50 iterations over
10 experiment realizations. Shades denote the standard deviation across
experiments.

Each type of agent utilizes one Deep Q-Network for maximizing rewards [15].
The rewards of the game, $R$, are as following:

  * •

Segregation reward. This incentive promotes agents to self-segregate. An agent
is rewarded +1 for each agent of similar kind that joins its observation
window, and -1 for each agent of different kind.

  * •

Interdependence reward. This incentive promotes interactions among agents of
different kinds. When an agent meets another agent of different kind, we
randomly choose a winner of the interaction (following hunting dynamics). The
winner (hunter) receives a positive reward, that we vary across experiments,
and an extension of its lifetime by one iteration.

  * •

Vigilance reward. This incentive promotes agents to stay alive by providing
+0.1 reward for every time step they survive.

  * •

Death reward. This incentive rewards negatively agents who die or are hunted
by agents of opposite kind. Agents receive -1 reward when they die.

  * •

Occlusion reward. This incentive rewards movements towards occupied cells
negatively. If an agent tries to move towards an occluded area, the agent
receives -1 reward.

  * •

Stillness reward. This incentive promotes the exploration of space. Agents who
choose to stay still receive -1 reward.

Every agent takes one action at each iteration. The sequence of agents who
take actions is chosen randomly. There are five possible actions for agents:
to stay still or to move left, right, up or down. Agents are confined to the
borders of the grid and cannot move towards agents of their own kind. If an
agent moves to a location occupied by an agent of the opposite kind, it
receives the interdependence reward and the opponent receives the death
reward.

Mathematically, agents of type A are represented as $-1$, B as $+1$, empty
space as $0$ and border as $-2$ on the grid. Hence every agent’s spatial
observation at time $t$ is $O_{spatial,t}^{(i)}\in
P^{11x11}\leavevmode\nobreak\ |\leavevmode\nobreak\ P\in\\{-2,-1,0,1\\}$.
Moreover, every agent has the information of its remaining normalized life
time, represented as $O_{age,t}^{(i)}\in R$. Full observation of the agent $i$
at time $t$ is $o_{t}^{i}\in O_{t}^{(i)}=O_{spatial,t}^{(i)}\cup
O_{age,t}^{(i)}$. Let $\phi_{A}$ and $\phi_{B}$ denote the Q-Networks of type
A and B. Then the networks’ goal is to satisfy Equations 1 and 2.

$\phi_{A}^{*}=\arg\max_{\phi_{A}}\mathbb{E}[\sum_{t=0}^{T}\sum_{i=1}^{N_{A}}\gamma^{t}r_{t}^{(i)}|o_{t}^{(i)}]$
(1)
$\phi_{B}^{*}=\arg\max_{\phi_{B}}\mathbb{E}[\sum_{t=0}^{T}\sum_{i=1}^{N_{B}}\gamma^{t}r_{t}^{(i)}|o_{t}^{(i)}]$
(2)

where $N_{X}$ denotes the number of agents of type $X$, $\gamma$ denotes the
discount factor, $r_{t}$ denotes the reward at time $t$ and $Q_{\phi_{X}}(.)$
denotes the Q-Network of agents of type $X$.

Each network is initialized with the same parameters. In order to homogenize
the networks’ inputs, we normalize the observation windows by the agents’ own
kind, such that positive and negative values respectively represent equal and
opposite kind for each agent. Actions are taken by following $\epsilon$-Greedy
exploration strategy. Exploration rate decays exponentially. In order to
stabilize the learning process, we use Adam optimizer [24], Experience Replay
[25] and Double Q-Learning [26]. Networks are trained in parallel over 12 CPUs
using data parallelism. We run one episode per experiment. Each episode is
comprised of 5000 iterations. Each experiment is repeated 10 times for
statistical analysis. Network details are given in Figure 1 (bottom) and
training details are given in Table 1.

##  3 Results

Figure 4: Dynamics of the experiment results. Panel (a) shows the average age
of agents. Panel (b) shows the maximum age of agents. Panel (c) shows the
percentage of agents that hunt at each iteration. Panel (d) shows the hunters’
cluster size prior to hunting the opposing agent. Panel (e) shows the age of
the agents hunting. Color is proportional to the interdependence reward (IR).
Darker color indicates higher interdependence reward. Each plot is obtained by
averaging results of ten experiments and 100 iterations. Shaded areas denote
the standard deviation among the experiments.

Experiments are conducted by setting up different values of incentives and
observing the emergent collective behavior associated with each experiment.
During simulations, agents explore the space of possible behaviors and inform
which behaviors are promoted under certain incentives and environmental rules.
As a result, we create an artificial environment for testing hypotheses and
obtaining information through simulations hard to anticipate given the
complexity of the space of possibilities.111Demonstration of the experiments:
(IR: 0) https://youtu.be/AgAeYMe2tUE (IR: 25) https://youtu.be/OZbl8qD50Mg
(IR: 50) https://youtu.be/Ca2p2cATmlw (IR: 75) https://youtu.be/R32Xu_EUpBQ.

In this case, we create agents who want to segregate from other kinds and
provide incentives to create interactions and interdependencies across kinds.
For this purpose, we model the Schelling dynamics for segregation and combine
it with the interdependence reward. The interdependence reward is given when
agents of different kinds compete and win against each other following hunting
dynamics. The one who is hunted dies and the hunter gets a positive reward and
life-extension. In total, there are four different experiments with
interdependence reward of 0, 25, 50 and 75 respectively. A set of videos are
available with one simulation for each setting. In the videos, colors
yellow/orange and cyan/magenta denote the types of agents. The color
brightness indicates the age of agents for both kinds.

Interdependence rewards diminish spatial segregation among different types. In
Figure 2a we show the collective behavior of the population, using heat maps
proportional to the probability of agents location during simulations
according to their type. The heat maps are visualized over one trial of the
experiments. Blue and red regions show biases towards each kind. White regions
show uniform occupation. The dynamics of segregation quickly result in patches
of segregated groups (top panels). As interdependence rewards increase, the
probability of one grid being occupied by agent of type A or B becomes uniform
and plots become white (bottom right panels). By creating interdependencies
among agents, they increase their interactions and reduce the spatial
segregation.

We measure segregation among agents using multiscale entropy. We convolve the
grid space with low pass filters of size 6x6, 12x12 and 25x25 using sliding
windows whose output is the window average value. We measure the entropy of
the distribution of window averages after each convolution across all
iterations. The segregation per iteration is defined as the average entropy
across the distributions resulting from the different filter sizes. The
resulting segregation dynamics is visualized in Figure 3. Segregation is high
when interdependencies are not rewarded (yellow curve). As interdependencies
increase (purple and black curves), the agents mix and the spatial segregation
is significantly reduced ($p<0.001$, see Section S2 in the Supplement).

Interdependencies affect the group dynamics. As we increase the reward for
interdependencies, the initially stable patches emerging from promoting
segregation become dynamic and mix with the other kind. The properties of the
population and associated activities reflect the change of dynamics. Agents
create an internal hierarchy where younger agents go out and hunt and elder
agents segregate and ensure reproduction. Evidences of such behavior are that
the average age of agents decreases and the hunting rate increases (Figure 4a
and 4c) and average hunter age is much lower than the average agent age
(Figure 4a and 4e). Moreover, the maximum age of agents per kind increases
(Figure 4b) showing that some agents stay protected and do not hunt. The
hunting strategy of agents is also affected by increasing interdependencies.
Pack size increases consistently with interdependence rewards. Figure 4d shows
the size of hunting clusters one step before hunting an agent. The increasing
cluster size given interdependence rewards suggests that agent association
yields better results. It also shows that hostile systems favor agglomeration
of agents for safety which can result in ultimate polarization. Additionally,
we also analyzed the effects of the vigilance rewards on the dynamics for
multiple reward values. Results show that higher vigilance reward increases
intra-kind interaction and results in more segregation (see Section S3).

Diverse areas attract younger people and people are older in segregated areas.
We show that older agents are more segregated than younger ones in the model
(see Figure 2b). The behavior has been observed in the model and verified with
human behavior using Census data. We analyzed the relationship between age and
segregation using Census data across the whole US (see Section S4). A
segregation metric based on racial entropy correlated positively with median
age by census tract (r=0.4). Our simulation shed light on an observation that
is not trivial about current societies.

In summary, our experiments show that increasing interdependencies among kinds
can be applied to reduce segregation. Moreover, hostile interdependencies will
result in in-group cooperation for hunting and competition for sheltering. The
emergent behavior of the population can be framed in the exploiter and
explorer discussion. A part of it chooses to segregate and another one to go
out and explore. The one who explores hunts and is vulnerable to be hunted,
but creates spatial integration. The one who segregates lives longer and
ensures reproduction of its own kind. In this model, explorers tend to be
younger and keepers tend to live longer. Spatial mixing was achieved by
increasing interaction rewards but was accompanied by larger clusters of
agents of the same size. Polarization may arise when there is an adversarial
relationship between the parts that segregate from each other. More generally,
emergent behaviors lie in a non-linear space where interaction properties
determine outcomes which may happen simultaneously and in different
combinations.

##  4 Discussion

We created an artificial environment for testing rules of interactions and
incentives by observing the behaviors that emerge when applied to multi-agent
populations. Incentives can generate surprising behaviors because of the
complexity of social systems. As problems become complex, evolutionary
computing is necessary to achieve sustainable solutions. We combine system
modeling (ABMs) with artificial intelligence (RL) in order to explore the
space of solutions associated to promoted incentives. RL provides ABMs the
information processing capabilities that enables the exploration of strategies
that satisfy the conditions imposed by the interaction rule. In turn, ABMs
provide RL with access to models of collective behavior that achieve emergence
and complexity. While ABMs provide access to the complexity of the problem
space, RL facilitates the exploration of the solution space. Our methodology
opens a new avenue for policy makers to design and test incentives in
artificial environments.

## Acknowledgements

We would like to thank Intel AI DevCloud Team for granting access to their
cloud with powerful parallel processing capabilities. Also, we would like to
thank Dhaval Adjodah for his valuable suggestions on training RL algorithms.

## Authors Contributions

ES, YBY and AJM contributed equally in the conceptualization, development and
interpretation of the experiments as well as in the paper write up.

## Data Availability

The source code of the model implementation as well as the data generated to
create this report will be made available upon publication.

## Additional Information

We declare that have no competing interests.

## References

  * [1] Lazer, D. _et al._ Computational social science.  _Science_ 323, 721–723 (2009).
  * [2] Eagle, N., Pentland, A. S. & Lazer, D.  Inferring friendship network structure by using mobile phone data.  _Proceedings of the national academy of sciences_ 106, 15274–15278 (2009).
  * [3] Morales, A. J., Vavilala, V., Benito, R. M. & Bar-Yam, Y.  Global patterns of synchronization in human communications.  _Journal of The Royal Society Interface_ 14, 20161048, DOI: 10.1098/rsif.2016.1048 (2017).
  * [4] Morales, A., Borondo, J., Losada, J. C. & Benito, R. M.  Measuring political polarization: Twitter shows the two sides of venezuela.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 25, 033114 (2015).
  * [5] Vosoughi, S., Roy, D. & Aral, S.  The spread of true and false news online.  _Science_ 359, 1146–1151 (2018).
  * [6] Ashby, W. R.  Requisite variety and its implications for the control of complex systems.  In _Facets of systems science_ , 405–417 (Springer, 1991).
  * [7] Taleb, N. N.  Black swans and the domains of statistics.  _The American Statistician_ 61, 198–200 (2007).
  * [8] Sayama, H.  _Introduction to the modeling and analysis of complex systems_ (Open SUNY Textbooks, 2015).
  * [9] Schelling, T. C.  Dynamic models of segregation.  _Journal of mathematical sociology_ 1, 143–186 (1971).
  * [10] Epstein, J. M. & Axtell, R.  _Growing artificial societies: social science from the bottom up_ (Brookings Institution Press, 1996).
  * [11] Axelrod, R.  A model of the emergence of new political actors.  In _Artificial Societies_ , 27–44 (Routledge, 2006).
  * [12] Cederman, L.-E.  _Emergent actors in world politics: how states and nations develop and dissolve_ , vol. 2 (Princeton University Press, 1997).
  * [13] Axelrod, R.  The dissemination of culture: A model with local convergence and global polarization.  _Journal of conflict resolution_ 41, 203–226 (1997).
  * [14] Epstein, J. M.  Agent-based computational models and generative social science.  _Complexity_ 4, 41–60 (1999).
  * [15] Mnih, V. _et al._ Human-level control through deep reinforcement learning.  _Nature_ 518, 529 (2015).
  * [16] Heess, N. _et al._ Emergence of locomotion behaviours in rich environments.  _arXiv preprint arXiv:1707.02286_ (2017).
  * [17] Sert, E., Sönmez, C., Baghaee, S. & Uysal-Biyikoglu, E.  Optimizing age of information on real-life tcp/ip connections through reinforcement learning.  In _2018 26th Signal Processing and Communications Applications Conference (SIU)_ , 1–4 (IEEE, 2018).
  * [18] Lanctot, M. _et al._ A unified game-theoretic approach to multiagent reinforcement learning.  In _Advances in Neural Information Processing Systems_ , 4190–4203 (2017).
  * [19] de Cote, E. M., Lazaric, A. & Restelli, M.  Learning to cooperate in multi-agent social dilemmas.  In _Proceedings of the fifth international joint conference on Autonomous agents and multiagent systems_ , 783–785 (ACM, 2006).
  * [20] Leibo, J. Z., Zambaldi, V., Lanctot, M., Marecki, J. & Graepel, T.  Multi-agent reinforcement learning in sequential social dilemmas.  In _Proceedings of the 16th Conference on Autonomous Agents and MultiAgent Systems_ , 464–473 (International Foundation for Autonomous Agents and Multiagent Systems, 2017).
  * [21] Sandholm, T. W. & Crites, R. H.  Multiagent reinforcement learning in the iterated prisoner’s dilemma.  _Biosystems_ 37, 147–166 (1996).
  * [22] Wunder, M., Littman, M. L. & Babes, M.  Classes of multiagent q-learning dynamics with epsilon-greedy exploration.  In _Proceedings of the 27th International Conference on Machine Learning (ICML-10)_ , 1167–1174 (Citeseer, 2010).
  * [23] Zawadzki, E., Lipson, A. & Leyton-Brown, K.  Empirically evaluating multiagent learning algorithms.  _arXiv preprint arXiv:1401.8074_ (2014).
  * [24] Kingma, D. P. & Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * [25] Lin, L.-J.  Self-improving reactive agents based on reinforcement learning, planning and teaching.  _Machine learning_ 8, 293–321 (1992).
  * [26] Van Hasselt, H., Guez, A. & Silver, D.  Deep reinforcement learning with double q-learning.  In _Thirtieth AAAI Conference on Artificial Intelligence_ (2016).
  * [27] Nikolov, N., Kirschner, J., Berkenkamp, F. & Krause, A.  Information-directed exploration for deep reinforcement learning.  _arXiv preprint arXiv:1812.07544_ (2018).
  * [28] Tang, H. _et al._ # exploration: A study of count-based exploration for deep reinforcement learning.  In _Advances in Neural Information Processing Systems_ , 2753–2762 (2017).
  * [29] Fu, J., Co-Reyes, J. & Levine, S.  Ex2: Exploration with exemplar models for deep reinforcement learning.  In _Advances in Neural Information Processing Systems_ , 2577–2587 (2017).
  * [30] Bansal, T., Pachocki, J., Sidor, S., Sutskever, I. & Mordatch, I.  Emergent complexity via multi-agent competition.  _arXiv preprint arXiv:1710.03748_ (2017).
  * [31] Macy, M. W. & Willer, R.  From factors to factors: computational sociology and agent-based modeling.  _Annual review of sociology_ 28, 143–166 (2002).

## Supplement

##  S1 Future Work

There are many potential improvements to our work. We classify directions of
future work under three categories: representation, training and
experimentation. Our method can be advanced by representing agents more
realistically such as introducing heterogeneous personalities to agents or
facilitating network structure over agents to promote alliances. Moreover,
training RL agents yield better results with sophisticated exploration
strategies [27, 28, 29]. In addition to exploration strategies, MARL is shown
to perform better with curriculum learning [30]. Our aim is to extend the work
on multi agent curriculum learning to our problem.

Schelling and Predator - Prey models cover just a little portion of the ABM
domain [31]. We are currently working on extending this artificial environment
to other ABMs, i.e. Axelrod model [13]. Our goal is to develop an easy
interface where policy makers and AI researchers can collaborate on solving
societal problems.

##  S2 Statistical Significance

We validate the significance of the patterns we observe along the execution of
the simulation as we change the IR incentive in Figure 3. We analyze the
distribution of values across the last 1000 interactions for each IR values
and test the difference among their averages. In Table S1 we summarize the
results of the statistical tests. The differences in averages are
statistically significant ($p<0.001$) across all pairs of curves.

##  S3 Vigilance Reward

We analyze the effects of the Vigilance Rewards (VR) on the dynamics of
agents. In Figures S1 and S2 we the impact on segregation and age distribution
for multiple values of VR. The results show that increased VR increases intra-
kind behavior and as a results increases segregation. Therefore, segregation
may also be fostered by other types of behaviors.

##  S4 Agents Age

We analyze the significance of the spatial distribution of agent ages shown in
Figure 2b. In Figure S3 we show the entropy of the spatial distribution of
agent ages at multiple iteration times and interdependence reward (IR),
together with a randomized case for comparison. The randomized case is
constructed by drawing agent ages on the grid from a uniform distribution
between (0, 1) and calculating the entropy. The difference between the random
case and the empirical results is significant in all cases. We tested
significance by comparing each curve. A summary of the test results are
presented in Table S2.

The behavior has been observed in the model and verified it with human
behavior using Census data. We analyzed the relationship between age and
segregation using Census data across the whole US. A segregation metric based
on racial entropy correlated positively with median age by census tract
(r=0.4). In Figure S4 we present a scatter plot of the segregation metric
(x-axis) and average age (y-axis) of each census tract (dots).

$R_{1}$ | $R_{2}$ | t-value | p-value
---|---|---|---
0 | 25 | 4.373 | 6.44e-06
0 | 50 | 15.935 | 0.0
0 | 75 | 23.890 | 0.0
25 | 0 | -4.373 | 6.44e-06
25 | 50 | 11.267 | 0.0
25 | 75 | 20.331 | 0.0
50 | 0 | -15.925 | 0.0
50 | 25 | -11.267 | 0.0
50 | 75 | 11.691 | 0.0
75 | 0 | -23.890 | 0.0
75 | 25 | -20.331 | 0.0
75 | 50 | -11.691 | 0.0
Table S1: Statistical tests results comparing segregation outcomes for different values of Interdependence Reward denoted as $R_{1}$ and $R_{2}$. We tested the difference in segregation over the last 1000 iterations for every pair of curves. The tests show that the average segregation differs across curves of different rewards. Figure S1:  Spatial distribution of agent types with varying Vigilance Reward (VR) (vertical) and Iteration (horizontal). Color indicates concentration of each type.  Figure S2: Spatial distribution of agent ages with varying Vigilance Reward (VR) (vertical) and Iteration (horizontal).  Figure S3: Entropy of Spatial distribution of agent ages with iteration. Random case is constructed by drawing agent ages on the grid from a uniform distribution and calculating the entropy. The difference in entropy corresponds to spatial clustering of ages. Entropy is normalized by the logarithm of the grid size. $V_{1}$ | $V_{2}$ | t-value | p-value
---|---|---|---
Random | IR: 75 | 18.755 | 0.0
Random | IR: 50 | 17.105 | 0.0
Random | IR: 25 | 16.347 | 0.0
Random | IR: 0 | 15.665 | 0.0
IR: 75 | Random | -18.755 | 0.0
IR: 75 | IR: 50 | 2.676 | 0.010
IR: 75 | IR: 25 | 3.772 | 0.0
IR: 75 | IR: 0 | 3.520 | 0.001
IR: 50 | Random | -17.105 | 0.0
IR: 50 | IR: 75 | -2.676 | 0.010
IR: 50 | IR: 25 | 1.214 | 0.231
IR: 50 | IR: 0 | 1.039 | 0.304
IR: 25 | Random | -16.347 | 0.0
IR: 25 | IR: 75 | -3.772 | 0.0
IR: 25 | IR: 50 | -1.214 | 0.231
IR: 25 | IR: 0 | -0.144 | 0.886
IR: 0 | Random | -15.665 | 0.0
IR: 0 | IR: 75 | -3.521 | 0.001
IR: 0 | IR: 50 | -1.039 | 0.304
IR: 0 | IR: 25 | 0.144 | 0.886
Table S2: Statistical significance of comparing the average age entropy for
different Interdependence Reward, denoted by $V_{1}$ and $V_{2}$. Random
curves are denoted by sampling entropy value per iteration from a uniform
distribution between zero and one. Figure S4: Age and racial segregation by
census tract (dots). Pearson correlation r annotated in the Figure.