Table Of ContentTransferable measurements of Heredity in models of the Origins of Life
Nicholas Guttenberg and Matthieu Laneuville
Earth-Life Science Institute, Tokyo Institute of Technology, Tokyo, Japan
Melissa Ilardo
Centre for GeoGenetics, Natural History Museum,
University of Copenhagen, Copenhagen, Denmark
Nathanael Aubert-Kato
Department of Information Science, Ochanomizu University, Tokyo,
5
Japan and Leading Graduate School Promotion Center, Ochanomizu University, Tokyo, Japan
1
0
We propose a metric which can be used to compute the amount of heritable variation enabled
2
by a given dynamical system. A distribution of selection pressures is used such that each pressure
n selects a particular fixed point via competitive exclusion in order to determine the corresponding
a
distribution of potential fixed points in the population dynamics. This metric accurately detects
J
the number of species present in artificially prepared test systems, and furthermore can correctly
8 determine the number of heritable sets in clustered transition matrix models in which there are
no clearly defined genomes. Finally, we apply our metric to the GARD model and show that it
]
accurately reproduces prior measurements of the model’s heritability.
E
P
.
o INTRODUCTION els which appear to give rise to forms of heredity[7–10].
i One of the difficulties in working with abstract models
b
- of something as broad as the Origins of Life, however,
Studies of the Origins of Life often come down to de-
q is that it is often unclear how to compare the results of
[ bate on the proper definition of ’life’. However, under-
thesedifferentmodelsbecausetheyuseverydifferentun-
lying that philosophical point are a number of distinct
1 derlyingmetaphors—chemicalreactionnetworks,math-
phenomena which are exhibited by modern biological
v ematical functions, combinatoric assemblages of materi-
5 life whose emergence must be explained. By studying
als, and cellular automata among others. In order to
4 these particular, quantifyable phenomena, it is possible
help connect those abstract efforts to the more concrete
7 tomakeprogressinourunderstandingeveniftheoverar-
1 goals of understanding the Origins of Life, we need de-
ching philosophical question remains murky. For exam-
0 vise a metric which can be applied very generally to a
ple, byfocusingon’replication’indetailtherehavebeen
. wide array of abstract models. Furthermore, we would
1 manyadvancesinunderstandingautocatalyticchemistry,
0 like this metric to be equally applicable to experimental
both experimentally and theoretically[1–6], even if those
5 systems, sothattheabstractmodelscanbeevaluatedto
systems may not fully qualify as life.
1 see whether or not they are truly predictive.
: OneofthemilestonesofatheoryoftheOriginsofLife
v Heritable variation in modern biological systems pro-
i is the understanding of how to bridge the gap between ceeds from a well-defined information-carying molecule
X
abioticconditionsandtheonsetofevolution. Alongwith
and well-separated individuals which comprise a popula-
r replication,themechanicsofDarwinianevolutionrequire
a tion. In this case one can conceive of a very clear way to
thatthesystembecapableofheritablevariationandalso
understand heredity by directly comparing the DNA of
be subject to selective pressures. Selective pressures are
an organism to its parent. This leads to understandings
readilyavailable—chemicalreactionratesareingeneral
of mutation, the process of fixation, and also a way to
very sensitive to environmental parameters — but it is
understand natural selection as a population-level effect.
more difficult to obtain the ability to inherit variations.
Inaskingquestionsabouttheemergenceoflife,however,
In order to understand how chemical systems become
it is likely that we will need to deal with intermediate
capableofevolutionwemustthenask: whatarethenec-
systemsinwhichthingsarenotsoclearlyseparated. We
essary and sufficient elements for a system to be capable
maynotknowwhattheinformation-carryingcomponent,
of heritable variation?
or the information may be distributed across all the de-
Because this question is very broad, it is well-suited grees of freedom of the system. In addition, if asking
to abstract models, in which the dynamics of the system questions about pre-cellular life, there will not necessar-
may be simplified in order to determine what elements ily be a way to distinguish individuals in the population.
are truly necessary. There are several examples of mod- Wemayalsonotknowthetimescaleofasinglegeneration
2
of replication or be working with systems in which the infinitessimally. Theothermechanismisthatifthereare
dynamics are more continuous and distinct generations a number of states that are local maxima in the replica-
cannot be defined. tionrate,aninfinitessimalchangeinthereplicationrates
In order to deal with this, we attempt to find a def- canalterwhichstateistheglobalmaximum. Thismeans
inition of heritable variability that does not depend on thatasmallchangeinselectionpressuremaygiveriseto
therebeingaknowninformation-carryingmolecule,well- a discontinuous large change in the long-term dominant
definedindividuals,oraknowntimescaleassociatedwith state. This second mechanism is what lets us count the
theturnoverofgenerations. Assuch,wewillnotreferdi- number of local maxima even if the dynamics only ever
rectly to the idea of replication with a particular fidelity, finds the global maximum. Since the two mechanisms
but instead try to look at the long-time behavior of dy- behave qualitatively differently under a distribution of
namical systems. Heredity is then related to the prop- selection pressures, we can focus on distinguishing the
erty of history-dependence. If a system’s state remains different types of variation, and by doing so the results
dependent on its history at long times, that enables it to should be relatively insensitive to the particular choice
retainandrememberinformation. Forexample,asystem of selection pressures (so long as the distribution is wide
with bistability is capable of remembering a single bit of enough to find all the local maxima).
information. Therefore, the central idea of our metric of heredity
What distinguishes heredity in particular from any is that what we are measuring has to do with the set of
form of memory is that it is memory which is ampli- ways in which a given dynamical system can respond to
fied extensively via replication. A system with heredity, a distribution of selective pressures. The distribution of
as opposed to one with just memory, can contain multi- selection pressures connects with all of the different pos-
ple copies of one or more pieces of information, and each sible variable degrees of freedom of the system, but due
of those pieces of information is constantly being copied. to the presence of replication the results of this are pro-
This distinction is what allows systems with heredity to jectedontothesubsetofthosevariabledegreesoffreedom
undergo evolution via population dynamics. The con- whicharealsoheritable. Essentially,wecanmeasurethe
sequence of this is that even if the dynamics of a piece ’genotype’ hidden behind the set of phenotypes we ob-
of the system would be bistable or multistable, if the serve, by seeing correlations between phenotypes which
various pieces of the system are amplified at different persist over a variety of different selection pressures (in-
ratesthenatlongtimesandforlargesystemsthefastest cluding pressures which explicitly differentiate between
growing state will dominate. In population genetics this those phenotypes).
is the phenomenon of competitive exclusion[11]. As a What remains is to devise a computational method to
result, the number of heritable states we detect will be evaluatehowmanydifferentgenotypesarepossible,given
dependentonhowlongwemeasurefor—inotherwords, a set of outcomes. In general, real systems will not only
because we do not a-priori know the appropriate choice havefluctuationswhichmayblurthedistinctionbetween
oftimescale, wecannotrecoverauniquemeasurementof different species, but also will not have a clear separa-
how much heritable variation is possible. Different in- tion between parts of the system which are responsible
vestigators with different systems would need to appro- for the mechanics of replication and parts of the system
priately choose that timescale for their experiments, and which are responsible for holding heritable information.
there may not be a good way to say what that timescale An autocatalytic chemical system, for example, will be
should in fact be. This makes it hard to compare results accompanied by the decomposition byproducts of all of
across models. the members of the autocatalytic cycle. Those decom-
Ifhoweveroneperformedtheaboveexperimentinava- positionbyproductsmaythemselvesvaryinrelativecon-
rietyofdifferentenvironments(selectionpressures),each centration depending on the particular kind of selection
environment would potentially pick out a different dom- pressure that is applied, even within the same ’species’
inant state at long times. This is still a somewhat arbi- of core autocatalytic cycle.
trary choice, as it is necessary to pick a particular distri- In general, we will not be able to guarantee that any
bution of environments tuned to the given system, such algorithmwefindwillperfectlyidentifythevariationbe-
that the perturbation is strong enough to re-order the tweenheritablestatesversusthevariationwithinasingle
competingstatesbutnotsostrongastocompletelyalter heritablestate. Thisproblemisnothingnew—inbioin-
their structure. However, there is an advantage to this formatics, it is often difficult to precisely define a bacte-
approach, in that it separates out two distinct mecha- rial species, and instead a clustering classification called
nisms for the system having a different steady-state con- ’organizational taxonomic unit’ must be used. What we
figuration. One mechanism is that as the environment candoistotrytocharacterizethebehaviorofthemetric
is perturbed, there is a corresponding perturbation to when given ambiguous cases and understand the failure
the location of the long-time stable state. That is, when modes. Assuch,wewillhavesomeabilitytodevisetests
a particular reaction rate is altered infinitessimally, the to recognize when the metric cannot be applied accu-
resting concentration of each compound is also altered rately and to provide bounds on the error.
3
First we will explain the algorithm for computing our had the best performance both in terms of robustness,
heredity metric. Then, we will apply this metric to an ability to detect different heritable states, and failure
artificial data set with a known number of species, dis- modes which are least confusing. Out of the techniques
tributed in different ways and with different mutation we tried, we found the best results using PCA combined
rates, in order to show that it can correctly determine with an analysis of the eigenvalue spectrum, and so we
the number of species present and also to show what focus on that.
happens when it breaks down. Next, we will examine The reasoning behind using PCA to detect the num-
another example system in which we construct a transi- ber of heritable states is that if there are N different
S
tion matrix describing a set of overlapping autocatalytic heritable states distributed in a space of much higher di-
networks,andshowthatthemetriccanaccuratelydetect mensionality,weexpectthattherewillnotbecollinearity
thenumberofmodules. Wewillthenapplythemetricto betweenanypairsofstatesinthecasewherethereisany
the GARD model, a non-trivial model system that has randomness involved in the location of the states in the
a distributed, ’compositional’ heredity which has been high-dimensional space. As such, a set of N points de-
S
extensively studied[8, 12–15]. fines a subspace of dimension N −1, which is what the
S
PCAisdetecting. Iffewercomponentsareretained,then
one heritable state must be explained as a superposition
METHODS of the others, but as long as the space of features is of
muchhigherdimensionthanthenumberofdifferenther-
Ouralgorithmfordetectingpotentiallyheritablestates itable states, this projection will generally have an error
isasfollows,takingasinputthedatamatrixPandgiving whosemagnitudecomesfromtheunderlyingdistribution
asoutputthenumberofheritablestatesofthedynamics of the heritable states of the system rather than coming
N . Each row of the data matrix is a single observation from the fluctuations. By using PCA and looking at the
S
of the dynamical system of interest at long times, driven eigenvalue spectrum, we attempt to detect that differ-
by a particular selection pressure (where each row corre- ence.
spondstoadifferentselectionpressure),andeachcolumn TheinputtothePCAisamatrixofdataP,forwhich
corresponds to a feature of the system — these can be each column is some observable property of the system
binary features or scalars which have been normalized se(features)andwhereeachrowisthefinalsystemstate
with respect to eachother to have similar variance. under a particular randomly sampled selection pressure.
These features are things such as whether a particular
chemical is present beyond a certain threshold concen-
Data: Matrix of observations P taken at
Result: Number of heritable states N tration, whether a certain gene is present or absent, etc
S
begin —theymustbechosenbytheinvestigatorappropriately
Subtract the mean from each column of P; tothemodelinquestion. Ingeneral,becauseoftheprop-
Add copies of data rows for regularization; ertiesofPCA,featuresshouldbenormalizedwithrespect
Principal Component Analysis of P → N
eigs toeach-otherandhaveazeromeanoverthedistribution
eigenvalues λ in descending order;
i of the data.
Set the index j =N ;
eigs
repeat ThePCAthenfindsasetofnewmutually-independent
Computeλ = 1(cid:80)2j λ (λ =0ifi>N ); features which are linear combinations of the given ones
> j i=j i i eigs
Set the threshold T =α(λ −λ )+λ ; such that the covariance between features is zero. These
1 > >
Set j to be the index of the first eigenvalue featuresaresortedwithrespecttotheamountofvariance
whose value is less than T; in the data associated with each feature — this infor-
until j not yet converged; mation is contained within the eigenvalues. If there are
The measured number of heritable states is certain directions which encode the differences between
N =j+1;
S heritable states, the variance associated with those di-
end
rections will grow as more features are added. This is
because that variance comes from the distance between
As part of the development of this metric, we tried a heritablestatesinthefeaturespace. Asoneincreasesthe
numberofdifferenttechniques: recursivefeatureelimina- dimensionalityofthatspace, distanceswithinthatspace
tion to find the subset of minimally correlated features, will grow monotonically, as each added dimension corre-
information-theoreticmeasuressuchassequenceentropy sponds to an additional strictly-positive contribution to
andmutualinformation,clusteringalgorithms(K-means the total distance.
andAgglomerativeClustering),anddimension-reduction At the same time, the per-eigenvalue variance that is
techniques(PCA).Therearemanymethodsthatallseem due to random fluctuations (e.g. uncorrelated with se-
feasible, and we do not have a strong argument to pick lective pressures) will remain roughly constant because
oneabovealltheothersonthegroundsoffirst-principles. as new random variables are added, new eigenvalues are
Instead, we present the method which out of our test set also added at the same rate. So as long as there are a
4
sufficient number of features, we expect that we should ideatoexamineaconvergenceplottoensurethattheal-
be able to use the eigenvalue spectrum to distinguish gorithm is behaving consistently for whatever choice of
those that come from the system being in globally dis- α ends up being used.
tinctstatesandthosethatcomefromnoise. Similarly,as We have made a simple Python script using
moredataisadded,weexpectthealgorithmtoconverge. Scikit-Learn[16] which implements this metric given
In general, assuming a sufficient number of features, the a file containing features and observations available
set of heritable states is a smaller-dimensional subspace at http://www.github.com/ModelingOriginsofLife/
than the set of possible random fluctuations. As such, Heredity so that other researchers can easily apply it
thedistributionofthedatawillconvergemorequicklyin to their own data sets.
directions corresponding to heritable states than it will
inthedirectionscorrespondingtofluctuations,andsowe
expect an increasing contrast as we add more data. Consistency checks
Once there are sufficient features and data points, we
can plot the eigenvalues sorted by rank. In cases where Our metric works best when the system has a small
there are multiple attractors for the system state, these number of heritable states compared to the number of
show up as a number of large eigenvalues. There is then features and independent measurements provided. For a
usuallya large drop, followedbya largenumberof small small number of heritable states, the gaps are generally
eigenvaluesassociatedwiththefluctuationsofthesystem verywell-defined,butasthenumberincreasesitbecomes
around these attractors. Some example eigenvalue plots harder to distinguish them from the noise. The algo-
are shown in Fig. 5. While the gap can often be located rithm strictly cannot detect more heritable states than
by eye, we need the algorithm to reliably detect it in an there are features. For this reason, as many potentially
automated way. For this purpose, we use an iterative relevant features as possible should be provided in the
procedure to localize the gap where there is a significant data. Because of this limit, systems with combinatori-
separationbetweenthesetofeigenvaluesbeforethebreak cally large genetic spaces must be broken up before they
and the set of eigenvalues after the break. can be analyzed with our metric.
In general, the eigenvalue spectrum from PCA on the This can be done by controlling the distribution of
expectedkindofdatahasaquickly-decayingpartoffixed selection pressures. Our metric only detects heritable
length (from heritable variation) followed by a slowly- states which are influenced by variation in the selection
decaying tail whose length depends on the number of pressure. If one constrains the selection pressure to a
features and samples. By doing a local average of the low number of dimensions, that will isolate a particular
eigenvalues below the cutoff, it is possible to detect the subset of heritable states. One can then use different se-
relevant height of the noise floor generated by mutation lection pressure distributions to find different subsets of
and the like in order to subtract it out. Roughly speak- heritable states in the same system. This is analogous
ing, this procedure finds the point at which the local to measuring the heritability of a single gene at a time,
derivative begins to be significantly steeper than the av- rather than looking at all possible genes at once.
erage slope from the origin to that point — essentially, When directly applying selection pressure to a par-
its a form of cliff-detection that normalizes with respect ticular subset of features, the variance in those features
toabackgroundmeanandmeanslope. Theparameterα will scale differently than the variance coming from un-
controls the sensitivity of the algorithm. If α is set to be derlying noise. This can create an artificial signal that
small, then in general the algorithm will detect finer dif- appears to correspond to a set of heritable states where
ferences between heritable states; however, when there really all that it is detecting is the distribution of exter-
are many heritable states, the error tends to be in the nal impulses. As such, it is best to exclude any features
form of large over-estimation of the number of states in which are being directly driven by the externally vary-
thesystem. Ifαismadelarger,thenthealgorithmtends ing selection pressure from the analysis. For example,
to saturate and cannot detect more than a certain num- in a chemical system with a number of compounds one
ber of heritable states, but errors are made in the form form of restricted selection pressure would be to add re-
of under-estimation. actions which decompose certain combinations of com-
In general, by examining the convergence pattern of pounds. The compounds which can be directly decom-
the algorithm for a few cases (as in Fig. 1), it is possible posedbythisreactionshouldnotbeincludedasfeatures
to find an optimum value for α which produces the most in the analysis.
stable results with as little data as possible. We recom- Another possible solution is to consider the properties
mend that α = 0.1 be taken as a good initial parameter of systems not in terms of a fixed count of the number
based on the results of our test cases. Sensitivity to the of heritable states, but in terms of scaling laws. If we
choice of α may be an indication that the algorithm is were to imagine applying our metric to DNA sequences
failingtodetectthedifferencebetweenheritablevariation oflengthLbeingreplicatedviathepolymerasechainre-
andthebaselinefluctuations,andsoitisgenerallyagood action,thenthenumberofheritablestatesisexpectedto
5
scaleasexp(L). IfLislarge,wecannotexpecttosample Species Detection
allpossibilities,butifLissmallthenitisstillpossibleto
doso. Assuch,wecouldimaginefixingL(orfilteringthe
In this test case, we artificially prepare a population
results according to L) and then performing the analysis
of distinct heritable states in order to see whether or not
for different L values. In the DNA system, we might be
themetricproperlydetectshowmanystatesarepresent.
able to see that the number of heritable states is scaling
Here, we do not try to see heritability emerge from an
in an open-ended way as we increase L (so the heredity
underlying dynamic but instead impose it directly in a
is unlimited), and that it is growing exponentially with
population genetics framework such that each heritable
L (so the heredity is combinatoric in nature). On the
state is the equivalent of a distinct species of organism.
other hand, if we were to do something similar in other
We implicitly allow the process of competitive exclusion
dynamical systems, we might see a characteristically dif-
proceedtocompletionineachcase,sothatthefinalstate
ferent scaling.
of the system for each independent run is a particular
There may also be problems when there is no heredity
species picked out of the set of possible species, plus mu-
at all in the system. When this happens, all directions
tation.
have the same apparent variance and so the behavior of
A system state is comprised of a number N of binary
methodstodetectthe’jump’betweensetsofeigenvalues f
features. Each species is a particular system state, with
often becomes ill-defined. This particular error can be
a fraction f of bases set to 1, and the rest set to zero.
detected by examining the convergence curve of the al- 1
When generating a novel end-state, we pick a random
gorithm periodically as one accumulates more data rows
speciesandthenapplymutationtoitsbinaryfeaturevec-
— if there is only one species, the number of detected
tor — each feature has a probability m of being flipped.
species will keep increasing as one adds data, until it is
We generate a population of N such states, which we
equal to roughly 30% of the number of features. P
then analyze using our metric.
It is also possible to use a form of regularization to
We discuss a number of particular configurations of
helppreventthiserrorinthefirstplace. Todoso,wear-
this process, in order to try to see how the algorithm
tificially add a small number of data rows with very low
responds to various possible problems that might arise
variance(essentiallyanartificial’species’),andthensub-
inrealdata. Thebasecasewewillconsiderhasf =0.5,
tractonefromthedetectednumberofheritablestatesat 1
m=0.05,N =1000,N =1000,andavariablenumber
the end of the analysis. In our examples, we do this by f P
ofspecies(sampleduniformly). Thisisarelativelygentle
takingarandomdatarowfromthesystemandmaking10
case, as the average distance between the binary vectors
copies of it. This guarantees some degree of consistency
for two different species is d ≈ 22, whereas the average
with the underlying statistics of the data (which might
distance between members of the same species is d ≈
be violated if one were to use a completely arbitrary ar-
10, and so the various heritable states should be well-
tificial species). We find that this does not significantly
clustered. The performance of the algorithm in this case
harm convergence elsewhere, but strongly helps in the
isshowninFig.2a. Inthisbasecase,thecorrectnumber
case where there is actually only one heritable state in
of species is detected to within one species error out to
the system.
N = 150 — this means that in many cases there may
S
only be a few individuals of those species in the system.
We may also want to measure systems in which fluc-
RESULTS
tuations are very large compared to the systematically
inherited variance, and so we want to see what happens
We present a number of test cases in order to eval- when the mutation rate m becomes large in order to see
uate the performance of our metric. We start with a how the algorithm fails. Our second configuration is the
straightforwardcaseinwhichwedefineapriorithestruc- same as the base configuration, but with a much higher
ture of the heritable states, and see if our metric can mutation rate: m = 0.2. This means that the average
detect the number of heritable states built into the sys- distance within a species is d ≈ 19 — close to the aver-
tem. Inthesecondcase,wemodelthedynamicsofmuta- age distance between a random pair of species. As the
tionandreplicationwithatransitionmatrixandevaluate number of species increases, we expect to have several of
the algorithm on the resulting attractors — this lets us therandomspeciesendupbeingclosertogetherthanthe
test whether the method of using a distribution of se- size of their mutational haloes. The results are shown in
lection pressures works for picking out heritable states. Fig. 2b. For small species counts, the algorithm is still
Thethirdtestcaseisthegradedautocatalysisreplication accurate, but at larger species counts there is a point
domain (GARD) model[8, 13] of compositional heredity at which the algorithm fails to detect all of the species
in replicating vesicles, which lets us test our metric on a present. AnincreasedN ,correspondingtotheavailabil-
P
non-trivial system whose heredity properties have been ity of more data samples, allows the algorithm to detect
studied elsewhere. thecorrectnumberofspecies. Asaresult,aconvergence
6
features corresponding to all combinations of the former
cases (e.g. m = 0.05,2 and f = 0.05,0.5 (Fig .2g). In
103 1
general,thealgorithmperformsreasonablywellforthese
cases, tending to make errors in the direction of under-
estimating the number of species.
The final case we look at is when the species them-
selves are not uniformly sampled. Instead, we sample
nt102 the various species with probability such that:
u
o
c
s
e
eci (cid:40)
p (1−1/N )i−1/N i≤N −1
ured s pi = (1−1/NSS)NS−1 S i=NSS
s
ea101
M
Thismeansthatthepopulationratiobetweenthemost
α=0.075 commonandleastcommonspeciesgrowsasNS−1. Even
α=0.15 with a low mutation rate (m = 0.05), this appears to
No regularization
make the problem significantly harder when the number
α=0.1
N =40 of species to distinguish grows. This makes sense, as at
100 S some point the rarest species only have a handful of ex-
0 1000 2000 3000 4000 5000 6000 7000 8000
N amplesinthepopulationandarecomparabletovariation
P
duetomutationalnoise. Asmightbeexpected,thealgo-
rithmdegradesintheformofunder-detectingthecorrect
FIG.1. Convergence of the measured number of heri-
number of species. One complication however is that as
tablestatesasthenumberofdatapointsisincreased.
This shows the effect of regularization and of the choice of α the (apparent) mutational floor grows, this can interfere
on the convergence pattern in this particular test case. The withthedetectionofspeciesthatwerecorrectlydetected
parameters for these data are f1 =0.05,m=0.2,NS =40. previously. This case is plotted in Fig. 2h.
When applying this algorithm to a novel system in
which the behavior is not already well-understood, it is
plot (Fig. 1) showing the measured number of species importanttorecognizepotentialsignsthatthealgorithm
calculated as a function of the number of data rows NP may not be detecting the correct number of heritable
may be useful in determining the trustworthiness of the states due to an insufficient number of features or data
metric for a given system. points. Tothisend,wevariedthenumberoffeaturesand
Anothercommoncomplicationisthatthedistribution data rows in this simple case and measured the number
of values of the binary features is not uniform. If for of species that can be included before the algorithm’s
example there is a particular chemical produced in only performance degrades. Specifically, we look at the point
oneoftheheritablestates,thenafeatureassociatedwith atwhichtheerrorinthemeasurednumberofspeciesex-
that chemical would be zero most of the time but rarely ceeds 50%. The results of this are plotted for an easy
would be one. This would not simply be an improbable test case (f = 0.5,m = 0.05, homogeneous population)
1
fluctuation,butinsteadsignifiesthepresenceofapartic- and a difficult test case (f = 0.05,m = 0.05, hetero-
1
ularheritablestate. Intermsofthistestcase, theconse- geneous population) in Fig. 3 (upper and lower panels
quence of these rare but significant feature values is that respectively). In the case of a homogeneous population,
theaveragedistancebetweendifferentspeciesisreduced. theabilityofthealgorithmtoresolvedistinctspeciesap-
We look at two situations here: one in which f = 0.05 pearstobelinearinthenumberofdatapoints(requiring
1
and m = 0.05 (Fig. 2c) and one in which f = 0.05 and roughly 4N data points in the easy case) and sublin-
1 S
m = 0.2 (Fig. 2d). This effect does somewhat destabi- ear in the number of features (though the algorithm can
lize the algorithm (especially in the high mutation case), never detect more species than there are features due to
but in general is less harmful than simply having a high the properties of PCA). When the population is hetero-
mutation rate, even though for f = 0.05 and m = 0.05 geneous, then the scaling with the number of data rows
1
the average distance between species is the same as the seemstobemoreseverelylimited,dependentonthefunc-
average distance between members of a species. tional form of the distribution of heritable states.
The next configurations we look at is when the binary
featuresarenothomogeneous. Welookatthreeexamples
of this: one in which half of the features have f = 0.5 Transition Matrix Model
1
and the other half have f =0.05 (Fig. 2e); one in which
1
halfofthefeatureshavem=0.05andtheotherhalfhave In the previous case, we artificially constructed a spe-
m=0.2(Fig.2f);andoneinwhichwehavefourtypesof cific distribution of species. In real cases, however, we
7
150
m=0.05,f =0.5, homogeneous
120 a) b) 1000 1
90 200
60
30 800
0 s
119250000 c) d) eatures 60020 460 80 100120140160180 112600# of specie
ecies 6300 # of f 400 80 mum
Sp 0 axi
of 150 2 M
# 2000 40
d 120 e) f)
e
ur 90
as 60 40 0
Me 30 1000 m=0.05,f1=0.05, heterogeneous 32
0 28
150 28
120 g) h) 800
90 24es
60 8 eci
3000 30 60 90 120 150 0 30 60 90 120 150 eatures 600 28 1260# of sp
NS # of f 400 0 24 12mum
xi
a
8 M
FIG.2. Behavior of the heredity metric on the species 20
200
detection test case. For these plots, perfect detection 16 4
corresponds to the line y = x. These cases use α = 0.1, 12
NP =1000, Nf =1000, and 10 extra duplicate sequences for 4 0 4 0 0
regularization. a) f = 0.5,m = 0.05, b) f = 0.5,m = 200 400 600 800 1000
1 1
# of data points
0.2, c) f = 0.05,m = 0.05, d) f = 0.05,m = 0.2, e)
1 1
f = 0.05,m = 0.05,0.2, f) f = 0.05,0.5,m = 0.05, g)
1 1
f1 = 0.05,0.5,m = 0.05,0.2, h) f1 = 0.05,m = 0.05, het- FIG. 3. Confidence regions for detecting different
erogeneous population numbers of heritable states. Upper: detectionconfidence
regions (50% error) for the easy case (m = 0.05, f = 0.5,
1
uniformdistributionofspecies). Inthiscasethealgorithmis
mostly data-limited, accurately detecting at most a number
expect to use a distribution of selection pressures in or-
of species roughly equal to 0.25N . Lower: detection confi-
der to probe the heritable states of a system. We extend P
dence regions for a harder case (m=0.05, f =0.05, hetero-
the idea of the former test case to take into account dy- 1
geneousdistributionofspecies). Inthisregimethealgorithm
namicalprocessesandtheresponsetoselectionbyusinga isfeature-limited—morefeaturesareneededtodetectmore
transition matrix style approach to a replicative process. species.
Inthismodel, theoverallsystemisdescribedbyadistri-
bution of random walkers over a set of nodes. There is a
transitionmatrixwhichgovernsthemovementoftheran- core of a well-defined species — e.g. as one increases
dom walkers between nodes. Additionally, replication is the number of potential nodes, the ratio between the av-
implemented by way of a per-node replication rate, such erage population of the nodes belonging to the species
that each iteration more walkers may leave a node than and the average population of the nodes not belonging
initially entered it. This can be thought of as writing a to the species diverges to infinity. On the other hand, if
transition matrix where the columns are not normalized mutation is stronger than replication then this ratio be-
to 1. comesindependentofthetotalnumberofnodes(e.g. the
This sort of model is similar to quasi-species in pop- populationfillsthespaceofpossibilitiessemi-uniformly).
ulation genetics. In such systems, there is an error- In our case, we would like to build a transition matrix
threshold[17] with respect to the balance between mu- such that its’ overall structure allows for the existence of
tation(linksfromanodetoothernodes)andreplication clustersofnodeswhichmutuallyexceedtheerrorthresh-
(links from a node to itself). When the selective ad- old. This is a closer analogy to autocatalytic chemical
vantage of a node is sufficiently strong compared to mu- systems, in which a given compound will not necessarily
tational diffusion, then that node ends up forming the directly replicate itself, but will instead proceed through
8
a number of intermediaries (which may have their own
side-reactions, or be part of multiple autocatalytic cy-
cles).
Suchtransitionmatrixbasedsystemsarelinear,which
means that the dynamics can be solved directly by com-
puting the eigenvalues of the transition matrix. At long
times, the state of the system will always be dominated
by the largest eigenvalue, which will tend to not be de-
generate unless there is a strong symmetry in the transi-
tion matrix (e.g. something like a block-diagonal struc-
ture). This is equivalent to competitive exclusion. As
such, even if we build a transition matrix with multi-
ple species which replicate with fidelity, at infinite time
we should only observe a single system-wide state corre-
sponding to the largest eigenvalue. This means that if
we wish to detect the cluster structure of the transition
matrix strictly from state data, we cannot do so unless
weapplyadistributionofselectionpressures(whichtakes
theformofpermutationstothetransitionmatrix)oroth-
erwisestochasticallydrivethesystemstronglyenoughto
overcome the gaps between its large eigenvalues.
To produce this type of topology, we use the following
procedure: FIG. 4. Transition matrix network structure. Example
structure of the transition matrix with 100 nodes, N = 20,
L
1. Randomly assign each node membership to a par- and N = 3 visualized as a network. The darkness of edges
C
ticular cluster k (out of N total clusters), and a is proportional to the probability of a transition along that
i C
replicationrater =r +δrη ,whereη isarandom edge. Eventhoughtherearemorebetween-clusterlinksthan
i 0 i i
variable distributed uniformly between [−1,1] within-cluster links, within-cluster links are a factor of 1000
morelikelytobefollowedthanbetween-clusterlinks,meaning
2. For each node, create N random links to other that overall a random walker tends to stay within its current
L
nodes. If the same pair of nodes is chosen multiple cluster.
times, sum the weights of the links
3. Foreachsuchlink,assignaweight: (1+γ)/2ifthe of half of the nodes by about 1% in order to create a
nodesbelongtothesamecluster,or(1−γ)/2ifthe distribution of selection pressures. In addition, we allow
nodes belong to different clusters. thereplicationratesateachnodetovaryrandomly(uni-
formly) by ±1% each iteration to represent the effects of
4. Normalizetheweightsoflinksleavingeachnodeso
a noisy environment.
that they sum to r
i Weusenetworkswith1000nodes,N =200,r =1.1,
L 0
This produces matrices with clusters that still have δr = 0.03, γ = 0.999, and a variable number of clusters.
weak inter-cluster connections. An example of such a We perform 250 iterations per run over 1000 runs and
network with 100 nodes, N = 3 clusters, and N = 20 use the resulting population vectors as our data points
C L
linkspernodeisshowninFig.4. Wegenerallyhaveeach for the heredity metric. The results are shown in Fig. 5.
nodeconnectto20%oftheothernodesinordertoensure For small numbers of clusters, our algorithm success-
that there are few nodes which do not link to any other fully detects the cluster count. However, as the num-
member of their cluster (which becomes a factor as the berofclustersgrows,between-clustertransitionsbecome
number of clusters grows). moreaccessibleduetotheincreasingnumberofbetween-
Given a transition matrix T of this form, we can at- cluster links compared to within-cluster links. This re-
tempt to apply our metric by simulating the dynamics sults in the gap in the eigenvalue plot shrinking, which
of this system. Rather than simulate individual random makes it harder to detect the correct number of clus-
walkers, we work with a vector encoding their popula- ters (it is somewhat equivalent to what happens with a
tion. At each iteration, we apply the transition matrix large mutation rate in the previous test case). When the
to this vector p =Tp and then normalize the vector number of clusters is larger than 32 or so, the gap seems
t+1 t
to 1 (this has no effect on the dynamics, but helps keep to disappear entirely (though there is clearly still some
the numerics stable). Each run, we start with a random kind of structure in the eigenvalue plot) — the result is
population vector and create a permuted transition ma- that our metric detects only the steep slope between the
trix by randomly reducing the on-site replication rates first few eigenvalues, and predicts a very small number
9
that qualify as sufficiently high fidelity to be of interest.
ThesestatesarethenclusteredusingK-meansclustering,
100
N =1 withthesilhouettemetricusedtodeterminetheoptimal
70 C
60 NC=4 number of clusters. As such, it is an example of an inde-
ed50 NC=8 pendentmeasurementofheritablestatesthatwecanuse
10-1 Detect423000 NNCC==1362 fmoretcroicmsphaoruisldonb.eTabhleedtiostoinpcetriaotne wisitthhoauttinexpprliincictipklneoowulr-
10
edge of the pattern of replication of individuals.
e 0
nvalu10-2 0 10203N0C40506070 givWene finoll[o8w],tuhseinimg pspleemciefinctaβtiomnaotfritchees GprAovRiDdedmobdyelthaes
e
g
Ei authors(asthedistributionofvaluesinthematrixisim-
portant to the degree of heredity observed). However,
rather than using the repeated growth and division of a
10-3
singleindividualvesicle,wesimulateapopulationofvesi-
cles under a distribution of selection pressures and then
only permit our metric to observe the system-averaged
10-4 compositions. The idea here is to test our ability to use
0 10 20 30 40 50
competitive exclusion to extract out the individual-level
Eigenvalue Rank
heredity from an analysis of the system that is not able
to distinguish the boundaries between individuals. We
FIG. 5. Eigenvalue plot corresponding to transition
use populations of size P = 1, P = 10, and P = 100,
matrices with different numbers of clusters. Whenthe
number of clusters is small, there is a clearly-defined gap be- and simulate for approximately 100 generations (direct
tweeneigenvaluescorrespondingtoclusteridentityandeigen- examination of the timeseries data shows that even for
values corresponding to the fluctuations. As the number of P = 1000 the system has reached steady by 30 genera-
clustersincreases,thegapclosesanditbecomesmoredifficult tions, andsothislengthoftimeshouldgenerallybelong
to detect the structure of the population. Inset: Number of
enoughtoconsistentlyreachsteadystateforsmallerpop-
clusters detected by the algorithm versus actual number of
ulations).
clusters.
Selectionpressurewasaddedbyspecifyingaparticular
target composition (cid:126)t and then multiplying the elements
of clusters compared with the actual amount. of the affinity matrix β according to:
β →(1+σt )β (1)
GARD ij i ij
where σ measures the strength of the selection pres-
The GARD model[8, 13] is a model which has been sure. We find that selection pressures in the range be-
observedtohaveaformofcompositionalhereditythatis tween σ = 5 and σ = 50 seem to be sufficient to change
significantly different than the type of heredity obtained the dominant state without destroying the structure of
via information-carrying polymers[18]. In GARD, the the underlying dynamics (whereas for σ = 500 the un-
fundamentalobjectsofthesystemarevesicleswhichcan derlying states are clearly distorted). We use σ = 50 in
be made from combinations of different types of lipids. our simulations, as the larger we can safely make σ the
The vesicles grow via a catalytic process — each lipid more easily we can sample rarer heritable states.
type has some affinity with each other lipid type (given We measure the mean number of heritable states for
by the random affinity matrix β), and thereby controls populations of size 1, 10, and 100 averaged over 24 dif-
the rate at which new lipids are added to the vesicle. ferent random β matrices. The measured results appear
When the vesicle reaches sufficient size N , it divides, to be consistent with eachother to within the measure-
max
randomly distributing its lipids to two child-vesicles of ment errorbars. We find that the number of heritable
size N /2. As such, there is replication built into the statesdetectedisconsistentwiththenumberspreviously
max
model, but that replication takes place (in general) at reported — on the order of about 3 distinct heritable
a very low fidelity. As such, one of the key questions states on average. Furthermore, we are able to detect
of GARD is whether or not high-fidelity states will nat- the effect observed by the authors that the number of
urally emerge from the dynamics as a result of a form heritable states tends to decrease with the ratio of the
of selection predicated on the robustness of a dynamical vesicle size and number of distinct lipid types. The re-
attractor. sults of our simulations are plotted in Fig. 6.
To detect these high-fidelity replicators, the authors In order to check our intuition as to what the PCA
look at the similarity of a vesicle to its parent and its metric is measuring and verify that it is detecting real
children[8]. This is used to filter out a subset of states distinctions, we look in detail at the dynamics of the
10
5.0
0.15
P=1
P=10
P=100
4.5 0.10
s
e state4.0 nent 0.05 Initial state
heritabl3.5 Compo 0.00
# of 2nd
ge −0.05
a
er3.0
v
A
−0.10
2.5
−0.15
−0.2 −0.1 0.0 0.1 0.2
2.0 1st Component
101 102 103
N
max
FIG.7. FlowdiagramofGARDdynamicsinheritable-
state space. This figure shows the average time-dependent
FIG.6. Hereditablestatesasafunctionofvesiclesize. behavior of the system state for a particular choice of the
In the GARD paper, the authors observed the trend that as β matrix and a distribution of selection pressures, projected
the vesicle size Nmax was increased, the number of heritable onto the first two PCA eigenvectors. The thickness of the
states decreased[8]. We show that our algorithm can detect streamlinesindicatesthemagnitudeofthevectorsintheun-
thistrendsuccessfullybyplottingtheaveragenumberofher- derlying averaged vector field. As the mean vector field is
itable states measured as a function of Nmax. The points are the result of averaging a number of random walks, this can
slightlyoffsethorizontallyforclarity,butallaremeasuredat be thought of as the strength of the local bias. The underly-
the same values of Nmax: 24,50,100,200,500 inggreyscaleintensityshowsthelogarithmoftheprobability
density for the system state being found at that location.
DISCUSSION
Oneofthedifficultiesofsynthesizingtheresultsofab-
stract modelling efforts is that often the details of the
GARD model for a specific choice of β for which our model are sufficiently idiosyncratic that it is difficult to
algorithm detected the existence of three species when directly compare the results of one model to another,
N = 50. We project the trajectory of the global sys- or to experiment. A model that works in terms of a
max
temstateforalargepopulation(P =1000)asafunction synthetic artificial chemistry may not easily make direct
oftimeontothethefirsttwoeigenvectorsofthePCAand predictions about steady-state chemical concentrations
then measure the mean direction of flow (averaged over in a real chemical system, because the details of chem-
the distribution of selection pressures) on the resulting ical energies and reaction rates will likely be different.
2D space. We find that individual trajectories go to iso- This means that although some insight can be gained,
lated points in the space, but as we vary the selection it is hard to say concretely what the one system shows
pressure an overall structure emerges in which there are us about the other. One way to overcome this limita-
clearly multiple distinct regions of the phase space — tion is to devise measurements which integrate out the
which region the system eventually ends up in depends detailsofthesystembutcapturefeaturesthataresome-
on the particular selection pressure. One region corre- what representation-invariant. One could, for example,
spondstoanextendedlinearsubspace, wherethedetails show that increasing the temperature in the synthetic
of the selection pressure position the system state at dif- system causes a certain statistical characteristic of the
ferentplacesalongtheline. Mostofthetimethewinning distribution of chemical compounds to change, and then
stateissomewhereonthisline. Theotherregionismore ask whether or not the same trend is observed experi-
broadly distributed but appears to have a second saddle mentally.
point and a corresponding pair of distinct peaks. This Much of the work being done to understand the emer-
structure is shown in Fig. 7. genceoflifefromprebioticchemistryfacesthisdifficulty.
