Heuristic optimization algorithms in the study of biological networks

Heuristic optimization algorithms in the study of biological networks Riu Rodr´ ıguez Sakamoto Heuristic optimization algorithms in the study of biological networks Riu Rodr´ ıguez Sakamoto Memoria presentada como parte de los requisitos para la obtenci´ on del t´ ıtulo de Doble Grado en F´ ısica e Ingenier´ ıa de Materiales por la Universidad de Sevilla. Tutorizada por Prof. Mar´ ıa del Carmen Lemos Fern´ andez Junio, 2019 Contents Abstract 1 1 Introduction and objectives 2 1.1 Heuristicalgorithms................................. 2 1.1.1 Approachable problems . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Algorithms ................................. 2 1.2 GeneRegulatoryModels .............................. 7 1.3 Objectives ...................................... 9 1.4 Usedtools ...................................... 9 2 Optimization across time 11 2.1 ModelingMethodology............................... 11 2.2 Motivation...................................... 11 2.3 Genetic Algorithm (GA) implementation . . . . . . . . . . . . . . . . . . . . . 12 2.4 Statistics: Receiver Operating Characteristic . . . . . . . . . . . . . . . . . . . 14 2.5 Results ........................................ 15 2.5.1 Network#1 ................................. 16 2.5.2 Network#2 ................................. 18 2.5.3 Network#3 ................................. 20 3 Optimization across space 23 3.1 Motivation...................................... 23 3.2 Implementation ................................... 23 CONTENTS ii 3.3 Encoding....................................... 27 3.4 Decoding....................................... 28 3.5 Logicgates...................................... 30 3.6 Representations ................................... 32 3.7 GAoperations.................................... 32 3.8 Results ........................................ 34 4 Conclusions 39 Bibliography 40 Abstract This work uses Genetic Algorithms (GA) to study in silico Gene Regulatory Models (GRM) and the relation between their structure and the temporal and spatial behavior. We first present a structure estimation algorithm based on the work by Ando and Iba [2] evaluating the networks temporal expression pattern, and in the second part of the thesis we apply a similar GA to find a GRM able to reproduce L. Wolperts French Flag model [28], evaluating its spatial pattern through simulation using CompuCell3D. 1 Introduction and objectives 1.1 Heuristic algorithms 1.1.1 Approachable problems Problems whose definition can’t allow analytical approaches and whose solution space is too broad to allow a exhaustive search can be found across many fields of science, where vast amounts of information are to be treated. Big Data has rightfully found its application in physics, biology, material science, among other domains in the last couple of decades due to technological advancements and an increasing availability of digitized data. Heuristic algorithms are used to trade accuracy, precision or completeness for speed to give an approximate solution to this kind of problems. 1.1.2 Algorithms The following three heuristic algorithms have all an inspiration in natural phenomena: ant colonies, metal annealing and genetics. A common term in heuristic algorithms is the fitness function, that is, the metric that describes how close a given solution is to achieving the desired goal. Depending on the problem, it can be a simple analytical function or as complex as a result of a simulation. Ant Colony Optimization Ant Colony Optimization is a family of algorithms inspired by the behavior of ants in a physical environment and aims to solve discrete optimization problems. Many ant species communicate between individuals and with the environment not visually but with the secretion and detection of chemicals called pheromones. Ants use these pheromone trails for example to mark the path to food sources. Other ants can then follow this path and reinforce the path leaving their own pheromones. This positive feedback characterizes the system, where the probability of an ant following a certain path increases with the number of ants that chose the path before. There are multiple implementations and variations of this family of algorithms, e.g. the Max-Min Ant System or the Ant Colony System. The first developed algorithm, called Ant 1. introduction and objectives 3 Algorithm 1 Ant System pseudocode (applied to the Travelling Salesman Problem). m←number of ants repeat for all k= 0 to mdo .Construct ant solutions i←random initial point for all construction step do j←next point to visit applying a random proportional rule i←j end for end for for all E(i, j)do .Update pheromone τij ←(1 −ρ)τij .Pheromone evaporation for all k= 0 to mdo ∆τk ij ←(1/CkE(i, j)∈Tk 0otherwise end for τij ←τij +Pm k=0 ∆τk ij end for until termination-condition Figure 1.1: In the case of the Travelling Salesman Problem (TSP), iand jrepresent the cities to visit. For each k-th ant from a total of mants, a route is built step by step. jis set by a random proportional rule (eq. (1.1)). τij is the pheromone present between points iand j.ρ∈(0,1] is the pheromone evaporation rate. E(i, j)represents the edge or path from ito j.Ckis the length of the tour Tkbuilt by the k-th ant, computed as the sum of lengths of the arcs belonging to Tk. Thus shorter tours get greater amounts of pheromone. System (Algorithm 1), consists of two phases which repeat until an optimal solution is found: 1. Ant’s solution construction: Each ant traverses the search space forming a path. 2. Pheromone update: This is done once all the ants have finished their paths, and the amount of pheromone deposited by each ant is a function of the paths fitness. The random proportional rule is a probabilistic action choice rule associated with the following probability of an ant kto go from point ito point j: pk ij =[τij]α[ηij]β Pl∈N k i[τij]α[ηij]βif j∈ Nk i(1.1) where ηij = 1/dij is a heuristic value available a priori defined as the inverse of the distance dij,αand βare parameters which determine the relative influence of the pheromone trail τij and the heuristic information ηij, and Nk iare the neighbor sites that ant kin site ihasn’t visited yet. For a more detailed explanation, along with variants of the Ant System algorithm, see the work by Dorigo et al. [7]. 1. introduction and objectives 4 Their applications are broad: Many NP-hard problems can be classified into one of these four categories: •Routing: Like the Travelling Salesman Problem presented in Algorithm 1, these problems consist of visiting a set of locations in an optimal order. AntNet is an Ant Colony Optimization algorithm designed to solve routing problems in telecommunications networks introduced by Di Caro and Dorigo [6]. •Assignment: This task consists of assigning a set of items Ito a number of resources R under some constraints, which can be seen as a mapping f:I → R. The objective function is a function of the mapping. •Scheduling: Allocating a limited amount of resources to tasks over time. •Subset: The solution to these kind of problems consists of a subset of the available components subject to constraints. E.g. the Weight Constrained Graph Tree Partition Problem, where an undirected graph with weighted nodes and arcs has to be split into different trees -group of connected nodesall with their weight between a minimum and a maximum. This NP-hard problem is found for example in the design of telecommunication networks and is generally not approximable [7]. Simulated Annealing This algorithm takes its name from the physical process of steal annealing, that is, the controlled cooling of the metal where initially at high temperatures the atoms in the molten state are randomly located, and later when it gradually solidifies, the steel finds an energetic minimum. Ideally, a slow enough annealing leads to the ground state. A rapid quenching might induce some irregularities that are trapped in, resulting in a solid with higher final energy (local minimum). The different final states of the metal correspond to the feasible solutions found by the algorithm. The energy corresponds in this analogy to the fitness function, and ground state to the optimal solution of the problem. A fast cooling is analogous to a local search, where the solution converges to a local optimum. The general operation of the simulated annealing algorithm is very similar to that of a local search: it starts at a random point xcin search space S. It evaluates xcwith the fitness function. Finds a new point xnin the neighborhood of xcand evaluates it too. If xnhas a better fitness than xc,xnreplaces xc. Else, xnreplaces xconly with a probability P(xn) = exp(−δ/T), where δis the difference in fitness of xcand xn. The general structure is described in Algorithm 2. 1. introduction and objectives 5 Algorithm 2 Simulated Annealing pseudocode (minimization problem). t←0 initialize T > 0 select random xc∈S evaluate xc repeat repeat select new point xn∈Snear xc if eval(xn)< eval(xc)then xc←xn else if U[0,1) <exp(−δ/T)then xc←xn end if until termination-condition T←g(T, t) t←t+ 1 until halting-condition Figure 1.2: xcand xnare points of the search space. If the fitness of xnis lower than xc(lower thus better, as this is a minimization problem), it replaces it and searches for a new xnin the vicinity of the newly assigned xc. If not, it replaces xcwith xnonly if a random value between 0and 1is lower than exp(−δ/T), where δis the difference in fitness of xcand xn. Genetic Algorithm Genetic Algorithms (GA) are based on the mechanics of genetics and natural evolution, and were firstly developed by Holland (1975) in his book Adaptation in Natural and Artificial Systems [13] as a method to both understand the natural adaptation processes in living organisms and to design a mathematical framework common to other fields like economics (’optimal planning’), artificial intelligence and psychology (’learning’). GAs search iteratively for a solution of an optimization problem. Each candidate solution (also called an individual) is a set of properties coded in a 1D array of data called chromosome. The algorithm starts with a random population of these candidate solutions and applies various operations to mutate and alter the individuals in the population each generation, until an individual with an acceptable fitness value has been found or until the maximum number of generations has been reached (Algorithm 3). The basic operators are: •Mutation: An individual is mutated in a random point of its chromosome, resulting in a slightly different individual with a similar fitness value to its original. This operator introduces some variability that can act as a local search in solution space. In the case of a simple 1D-array of bits, a random position mis chosen and it’s value is flipped. •Crossover: Two individuals exchange parts of their chromosomes. The specific way of crossing two chromosomes depends on the problem design, but in the simple case of a 1D-array of bits, a random position pis chosen and each sub-array is exchanged with 2. optimization across time 12 0 0.2 0.4 0.6 0.8 1 Expression level 0 20 40 60 80 100 Population distribution (%) X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20 X21 X22 X23 X24 X25 X26 X27 X28 X29 X30 Nodes X1 X26 Figure 2.1: Histogram of stable expressions levels from random initial conditions. 2.3 Genetic Algorithm (GA) implementation We want to obtain the values of the weight matrix wji that give rise to a certain time expression pattern. To format the matrix in a way that a GA can use, we concatenate the rows of the matrix in a single 1D array of real values, as done in [2]. As we are using a GA, we have to define how we implement the fitness,mutation,crossover and selection operators explained in Section 1.1.2. Fitness The fitness of an individual is computed from two values: the linear norm δand the parsimony factor P. •The linear norm between the generated expression pattern xi(t)and the target expression pattern yi(t)through time is given as: δ= N X i=0 |yi(t)−xi(t)|(2.3) 2. optimization across time 13 •The parsimony factor Pis proportional to the number of pathways through the GRN and the number of non-zero elements mof the influence matrix. P=T·m(2.4) T=0.01 ·#pathways N2(2.5) The value 0.01 is simply a scaling factor to be combined later with the linear norm into a single fitness value (2.8). It is the same as the one used by Ando and Iba [2]. The number of pathways is strictly speaking infinite, because the GRN is a cyclic graph. We limit the number of pathways calculating them using the adjacency matrix. The adjacency matrix is a square matrix whose element Aij is 1if there’s an edge from node ito node j, and 0otherwise. Eis the set of all edges of the graph. Aij =(1if (i→j)∈ E 0if (i→j)/∈ E (2.6) If we raise Ato the n-th power, each element (An)ij has the number of paths from node ito node jconsisting on nedges or hops. Therefore, if we sum all elements of Anfor n up to the number of nodes N, we will cover all possible paths through the graph. #pathways = N X nX i,j (An)ij (2.7) The justification behind the use of this second metric Pis based on the Minimum Description Length principle [2]. Using P, the search is biased towards smaller, less interconnected networks, e.g. to the parsimonious model. The total fitness is a weighted combination of the two values. DEAP allows to give multiple fitnesses individually and their corresponding weights. fitness =w1·δ+w2·P(2.8) Two-point crossover The crossover is the standard two-point crossover: A randomly selected position splits each individual in two. One section of one individual is swapped with the same section of the other individual, maintaining the total length of both individuals. Iterating in pairs over the complete population, the probability of mating two individuals is pcx. Gaussian mutation Mutation is done with a set probability pmut for each individual of the population in a given generation. For an individual to be mutated, each of it’s elements has a probability pind to suffer a Gaussian mutation of mean µ= 0 and standard deviation σ= 0,2. The values set for these hyperparameters is justified, as it often happens in many optimization techniques, through trial and error. We also take these values from the ones used by Knabe et al. [16]. 2. optimization across time 14 Tournament selection Tournament selection is also a standard procedure: it selects the best individual among three randomly chosen ones and repeats until it obtains enough to make a new generation. The criterion to choose the best individual is their fitness. 2.4 Statistics: Receiver Operating Characteristic Strictly speaking, a receiver operating characteristic curve is used to evaluate a binary classifier. It plots the true positive rate against the false positive rate. These statistical measures can be defined in our case for comparing graphs following the criteria shown in Table 2.1. •True positive rate (TPR) is the proportion of actual activating edges in the target graph that are correctly identified as such in the acquired graph. It is also called sensitivity. TPR=TP P=TP TP +FN = 1 −FNP =sensitivity (2.9) where TP are true positives, Pare total positives, FN are false negatives and FNP is false negative rate. •False positive rate (FPR) is the ratio between the number of repressing edges wrongly categorized as activating and the total number actual repressing edges. It can also be defined as 1−specificity, where the specificity or false positive rate is the proportion of actual repressing edges that are correctly identified as such. TNR =TN N=TN TN +FP = 1 −FPR = 1 −specificity (2.10) where TN are true negatives, Nare total negatives, FP are false positives and FPRis false positive rate. Measure Definition TP An edge that is activating in both the target and acquired graphs. TN An edge that is repressing in both the target and acquired graphs. FP An edge that is repressing in the target graph but activating in the acquired graph. FN An edge that is activating in the target graph but repressing in the acquired graph. PTotal number of activating edges in the target graph. NTotal number of repressing edges in the target graph. Table 2.1: Statistical measures defined for graph comparison. 2. optimization across time 15 2.5 Results The Genetic Algorithm parameters used are different for each graph. In general, bigger or more complex GRNs require more individuals per generation (the population pop is higher) and the GA should run more generations to find a suitable solution. The parameters for each network are shown in Table 2.2. For the first target graph, the target and acquired networks are shown side by side in Figure 2.2. The evolution of the sensitivity and specificity (Section 2.4) is shown in the Receiver Operating Characteristic in Figure 2.3, and the fitness in Figure 2.4. The comparison of the expression levels evolutions across time is displayed in Figure 2.5. Similar results were obtained for network #2 (Figures 2.6, 2.7, 2.8 and 2.9). However, for bigger GRNs (network #3: Figures 2.10, 2.11 and 2.12), the algorithm fills the weight matrix ωij, even against the parsimony factor Pintroduced in the fitness function, resulting in overregulated GRNs. The results are summarized in Table 2.3. For network #1, the Genetic Algorithm successfully replicated 8 out of the 10 edges in the network, out of the 72= 49 possible, with a edge weight error of 14.96%. This edge weight error is determined from the average difference: e=1 n n X i=1 ω(t) i−ω(f) i(2.11) where n=E(t)is the number of edges in the target network, i.e. the sets cardinality, and ω(t) iand ω(f) iare the weights of the i-th edge in the target and final graphs respectively. If the final graph doesn’t have the corresponding i-th edge, ω(f) i= 0. #N pop gen pcx pmut pind 1 7 1000 150 0.99 0.01 0.01 2 10 7000 150 0.99 0.01 0.01 3 30 8000 600 0.99 0.01 0.01 Table 2.2: GA parameters for the tested networks. Nis the number of nodes, pop the number of individuals per generation in the GA, gen is the total number of generations to be computed. pcx is the probability of mating two individuals, pmut the probability of mutating an individual, and pind the independent probability of each attribute within an individual to be mutated. # Fitness Edges in target Edges in final Missing edges Extra edges e 1 0.012434 10 9 X1→X1X6→X6 14.96% X5→X6 2 0.006762 13 13 X8→X7X7→X1 4.35% 3 0.207424 35 73 see Fig. 2.12 and Fig. 2.11 31.36% Table 2.3: Properties of the final networks obtained. eis the edge weight error (2.11). 2. optimization across time 16 2.5.1 Network #1 X1 0.5 X3 0.5 X2 -0.8 X4 0.7 X5 0.5 -0.5 X6 0.1 X7 0.9 -0.8 0.4 (a) Target graph X1 X3 0.601 X2 -0.786 X4 0.735 X5 0.309 -0.468 X7 0.958 X6 -0.453 0.146 0.283 (b) Acquired graph Figure 2.2: Comparison between the target and acquired GRNs for network #1. Green edges indicate activating regulations (ωij >0) and red edges indicate inhibitory regulations (ωij <0). X1→X1and X5→X6are missing regulations, and X6→X6is an extra regulation. Figure 2.3: Receiver Operating Characteristic evolution for network #1. Center is random, top left is better. The initial population of GRNs starts at the center, and gradually moves towards better individuals as generations evolve. 2. optimization across time 17 Figure 2.4: Evolution of fitness (objetive function) across generations in scalar (left) and log (right) scales for network #1. Although the best individual (lowest curve, orange line) keeps scoring better fitness values until the end, it’s drop rate diminishes significantly, together with the average fitness (blue line) at around generation 50, which suggests a decrease in the marginal efficiency of the algorithm for this network. Figure 2.5: Comparison between the target (solid line) and acquired (dotted line) GRNs expression levels across time with identical initial conditions for network #1. A decent overlap of the curves can be observed, which supports the low fitness value obtained for the best individual. 2. optimization across time 18 2.5.2 Network #2 X1 0.6 X5 -0.4 X2 X3 0.4 X4 0.6 X6 -0.5 X7 0.4 -0.7 X8 0.5 -0.3 X9 0.3 X10 -0.7 -0.5 0.3 -0.1 (a) Target graph X1 0.487 X5 -0.344 X2 X3 0.425 X4 0.644 X6 -0.472 X7 0.113 0.386 -0.824 X8 0.498 -0.325 X9 0.302 X10 -0.68 -0.536 0.28 (b) Acquired graph Figure 2.6: Comparison between the target and acquired GRNs for network #2. Figure 2.7: Receiver Operating Characteristic evolution for network #2. Center is random (starting point), top left is better. 2. optimization across time 19 Figure 2.8: Evolution of fitness (objective function) across generations in scalar(left) and log (right) scales for network #2. Figure 2.9: Comparison between the target (solid line) and acquired (dotted line) GRNs expression levels across time with identical initial conditions for network #2. 2. optimization across time 20 2.5.3 Network #3 (a) Receiver Operating Characteristic evolution for network #3. (b) Evolution of fitness (objective function) across generations for network #3. Figure 2.10: ROE and fitness evolution across generations for network #3. 2. optimization across time 21 X1 X5 1.0 X6 1.0 X2 X7 0.5 X3 0.4 X4 X8 0.2 X11 0.4 X9 1.0 -0.1 X10 0.3 -0.2 X13 -0.6 X14 1.0 X15 0.2 X16 0.5 X12 -0.2 X17 0.5 X19 0.1 X20 0.7 X24 -0.2 X21 -0.2 X22 0.5 X23 0.2 X18 -0.10.3 X25 0.4 X26 0.6 0.4 0.1 X27 0.6 0.3 X28 0.5 X29 0.4 X30 0.6 0.1-0.2 Figure 2.11: Target graph representation of the GRN for network #3. 3. optimization across space 28 Binary 0000 0001 0010 0011 0100 0101 0110 0111 kaor kd0.0625 0.125 0.1875 0.25 0.3125 0.375 0.4375 0.5 n0 1 2 3 4 5 6 7 Binary 1000 1001 1010 1011 1100 1101 1110 1111 kaor kd0.5625 0.625 0.6875 0.75 0.8125 0.875 0.9375 1.0 n8 9 10 11 12 13 14 15 Table 3.2: Binary encoding of all possible values for ka,kdand n. 3.4 Decoding We decode the rate equations from the genome. In order to do this, we interpret that when a protein Pjactivates another Pi, this activation is a reaction catalysis in biochemistry: the activator protein Pjbinds to a specific location Sof the DNA and forms a complex Pi, allowing the transcription of Pi. The location of the DNA can either be free (Pjnot present) or bound (with Pj), so the total is Ptot =S+Pi. The reaction is: S+nPj ka  kd Pi(3.1) where Sis the substrate, Pjis the catalyst, Piis the product, kais the association constant and kdthe dissociation constant. Sis the molecule that forms the complex with nmolecules of Pj to result in Pi. The three constants ka,kdand nare the ones present in each CRM in Figure 3.2. As it is derived by U. Alon [1], the association rate is proportional to the concentration of Sand the concentration of Pjto the power n(the probability of finding nmolecules of Pj simultaneously). The proportionality factor is the association constant ka: association rate =kaSPn j(3.2) The dissociation rate is proportional by a factor kdto the concentration of Pi: dissociation rate =kdPi(3.3) In equilibrium, the rate of change of the concentration of Pi, that is, the difference between association and dissociation rates, is zero in a steady state approximation: dPi dt =kaSPn j−kdPi= 0 (3.4) From here, we can obtain the apparent dissociation constant defined by the Law of mass action as: (Kij)n=kd ka =S·Pn j Pi (3.5) Solving for the ratio between bound Piand total S+Piwe get the Hill equation: Pi S+Pi =Pn j Kn ij +Pn j (3.6) 3. optimization across space 29 In non-equilibrium, the transcription rate for a protein Pifollows the Hill kinetic equation, expressed as: dPi dt =vmax Pn j Kn ij +Pn j =vmax (Pj/Kij)n 1+(Pj/Kij)n(3.7) where nis the Hill coefficient and Kij is the apparent dissociation constant (3.5). The Hill kinetic equation (3.7) is also known simply as Hill equation, which is how it is going to be named from here on, but should not be confused with the similar (3.6), which is not a differential equation. For n= 1, (3.7) reduces to Michaelis-Menten kinetics, one of the best known models of enzyme kinetics: dPi dt =vmax Pj Kij +Pj (3.8) Analogously, the Hill equation for gene Pjrepressing Pihas the form: dPi dt =vmax 1−Pn j Kn ij +Pn j=vmax 1 1+(Pj/Kij)n(3.9) The general form of the rate equations is a combination of both and is called the generalized Hill equation, as proposed by Del Vecchio and Murray [5]: dPi dt =vmaxA({Pj}) + vb 1 + A({Pj}) + I({Pj})−βPiPi(3.10) where βPiis the decay rate of the i-th protein, vbis the basal transcription rate and vmax is the maximum transcription rate (saturation). The functions A({Pj})and I({Pj})are the sums of normalized concentrations of the activating and inhibitory proteins of Pito the power of their respective Hill coefficients ni, described by the CRMs of Pi. For example, if P1and P2activate P3, the activator function has the form: A(P1, P2) = P1 K31 n1 +P2 K32 n2 (3.11) We can assert that the transcription rate of Pidoesn’t get over vmax and that the sigmoidal shape of the Hill equation is conserved for limiting cases (A→0,A→ ∞,I→0,I→ ∞) (Figure 3.3). We can also assert that, by means of the limiting factor of the decay rate, the model doesn’t explode: The concentration of Pireaches an equilibrium at a finite positive value: Peq i=1 βPi vmaxA({Pj}) + vb 1 + A({Pj}) + I({Pj})(3.12) 3. optimization across space 30 10 5 10 4 10 3 0 Activator (A) 10 5 0.5 10 2 10 4 1 Inhibitor (I) 10 3 1.5 10 1 10 2 2 10 1 2.5 3 Figure 3.3: Generalized Hill equation without the decay factor (βPi= 0) as a function of activator (A) and inhibitor (I) concentrations, with parameters νmax = 3,νb= 1,KiA = 50, KiI = 50,nA= 1,nI= 1. For both low Aand I, the transcription rate dPi/dt is the basal rate νb. For high Aand low I, it saturates at νmax and for low Aand high Ithe transcription rate tends to zero. 3.5 Logic gates As illustrated in Figure 3.2, each CRM can have more than one transcription factor (TF) than can act in conjunction as inhibitor or activator of an influenced protein PR. This allows to model logic gates, primarily AND and OR gates. The implementation of logic functions in DNA has been experimentally shown by Okamoto et al. [21]. AND, OR, YES and NOT are the basic logic gates from which all logic functions can be constructed from [9]. We have designed a genome structure (Fig. 3.2) and a set of kinetic equations (3.10) flexible enough that certain combinations give rise to logic functions as an emergent property of the system, rather than hard-coding them explicitly for example as a special ’AND’ gene. Remember that genes represent regulations or edges between nodes -the proteins-. The resulting kinetic equations for logic functions are equivalent to those derived by thermostatistical modeling from principles of statistical physics [3]. We can always design a genome to have such behavior (Fig. 3.4): •AND: Modeled as multiple TFs per CRM. The rate equation will depend on the product 3. optimization across space 31 of the concentrations, creating an effective AND gate (extended to continuous values). dP3 dt = vmax P1 K13 nP2 K23 n +vb 1 + P1 K13 nP2 K23 n−β3P3(3.13) •OR: Modeled as multiple CRMs activating a PR in a gene. If any of the CRM activates, the formation of the corresponding PR will increase. dP3 dt = vmax hP1 K13 n +P2 K23 ni+vb 1 + hP1 K13 n +P2 K23 ni−β3P3(3.14) •YES: Modeled as a PR with a single activating CRM. dP2 dt = vmax P1 K12 n +vb 1 + P1 K12 n−β2P2(3.15) •NOT: Modeled as a PR with a single inhibitory CRM. dP2 dt =vb 1 + P1 K12 n−β2P2(3.16) AND genome 1 β 2 β 3 β 3 2 1 a kd k n pr1 pr2pr3 1 CRM 1 G OR genome 1 β 2 β 3 β 3 2 1 a kd k n pr1pr3 1 CRM 2 1 a kd k n pr2 2 CRM 1 G YES genome 1 β 2 β 3 2 1 a kd k n pr1pr2 1 CRM 1 G NOT genome 1 β 2 β 3 2 0 a kd k n pr1pr2 1 CRM 1 G AND 1 2 3 1 2 3 1 2 1 2 Figure 3.4: Minimal examples of GRNs implementing each logic gate, together with their corresponding graph representation. 3. optimization across space 32 3.6 Representations We have presented three spaces in which an individual or solution can be represented: •Genome space: Each point in this space represents a 1D-array genome. The genetic algorithm applies operations of crossover, mutation and selection in this space. •Graph space: Points in graph space represent the decoded GRNs. •Space of systems of differential equations: Each point in this space represents a system of differential equations governing the dynamics of the protein concentrations and interactions. There is a bijective correspondence between elements in graph space and space of systems of differential equations. However, the decoding function T(3.4) is not injective nor surjective, because it can transform many genomes to the same GRN - e.g, two genomes consisting of the same genes but in different order -, and because some GRNs can’t be encoded in a genome - e.g, a GRN where one protein has a higher decay rate than the bounded (0,1] range that a genome can express-. The genome space contains the full hereditary information of individuals, so it can be called genotype. On the other hand, the space of systems of differential equations and the graph space, together with environmental information like initial conditions for the simulation, comprise the individuals behavior and observable characteristics, so they represent the phenotype. 3.7 GA operations Mutation and crossover are designed so that every possible outcome is a valid individual. Therefore we avoid having to reject invalid individuals or defining a penalty function. Mutation Mutation is simple: It looks for a random position in the genome where there’s a 0or a 1, and flips it’s value. We avoid changing the value of positions where there’s a 2or a 3because these digits mark the structure of the genome (starting positions of a gene and a CRM respectively, as explained in 3.3), so modifying these values would result in a drastic change in phenotype. Crossover To mate two genomes, we first choose a random section in each genome. If the chosen section is the first (where the decay rates β1, . . . , βNare encoded), we choose a random position within the section and exchange each part of the section between the two individuals. If the chosen section is a gene (second section onwards), we choose a random position within the section 3. optimization across space 33 but also an offset (from a normal distribution with its width or standard deviation being the address bit length, 3) which we apply to the left of one individual and to the right of the other to locate the final cutting points. We then split the genes in each individual and exchange them. This is similar to the crossover mechanism proposed by Knabe et al. [16]. Selection We combine two selection algorithms as done in [16]: Tournament selection, where a fixed number of individuals (30 in this case) is randomly chosen from the population and the individual with the highest fitness score is selected, and elitism, where the best individual of the population is always conserved. X X X-d X+d Mutation P C 01001011301021001001110001001010 01001111301021001001110001001010 Crossover (first section) P1 P2 01001011301021001001110001001010 01100101300120011110010011011 C1 C2 01001101301021001001110001001010 01100011300120011110010011011 Crossover (gene) P1 P2 01001011301021001001110001001010 01100101300120011110010011011 C1 C2 01001011301021000011011 01100101300120011110011001110001001010 Figure 3.5: Mutation and crossover mechanisms introduce generic variation. Mutation requires one parent (P) and produces one child (C) varying one of its bits. Crossover takes two parents (P1 and P2) and produces two children (C1 and C2). First-section crossover exchanges the bits on the right of a random position from the first section of both genomes. Gene crossover swaps bits to the left of X−dof P1 and the bits to the right of X+dof P2, being dan offset chosen from a normal distribution3with width σ= 3. pop gen pcx pmut tournsize 500 150 0.90 0.01 30 Table 3.3: GA parameters for GRN optimization across space. pop is the number of individuals per generation in the GA, gen is the total number of generations to be computed. pcx is the probability of mating two individuals, pmut the probability of mutating an individual. tournsize is the number of individuals from which the best is selected in tournament selection. 3The standard deviation σis chosen to be the same as the protein address length. Again, for greater number of proteins in the GRN, this address length is increased so σwill also be greater. 3. optimization across space 34 3.8 Results Figure 3.6 shows the evolution of the evaluated function -fitnessand the genome length across generations. Each generation -along the abscissahas a population for which there’s an individual with a minimum fitness - the best individual, in orangeand an individual with the maximum fitness value - worst individual, in green-. We also calculate the average fitness of the population - blue line-. Similarly for genome lengths, we plot the minimum, maximum and average genome lengths across generations. Figure 3.6a shows that from about generation 120, the best score stabilizes at fitness ≃4.4. Figure 3.7 shows the final stable expression profile for the best individual after 150 generations along the x-axis of the 2D plane, averaged from the values in the y-axis. The green curve is the protein sensing the morphogen gradient pr12.pr00 and pr01 are the color determinants. It shows resemblance to the objective French flag model (fitness = 4.401). To reconstruct the French flag model, we have compared an ad hoc colormapping based on the spatial fields of pr00 and pr01 with the theoretical french flag in Figure 3.8. Finally, figure 3.9 shows the complete regulatory network connecting the eight nodes of the system. Proteins pr02 and pr03 don’t have any inputs, which is in accordance with the constant zerovalued field in Figure 3.10. As the fitness does not include any parsimony factor, there are two negative consequences: •Genomes tend to become longer each generation, adding complexity. •The final GRN contains some residual elements which do not contribute to the final expression pattern, like pr02 regulating an AND gate, effectively suppressing its output as the output of an AND gate depends on the product of input concentrations and pr02 has always zero expression level. However, this information sparsity is desirable because of two reasons: •Unused sections of the genome can be put into use once again in the next generation after population crossover or mutation. •It adds diversity to the individuals of the population, observable in the wide range of fitnesses that the evolution encompasses - difference between maximum and minimum fitnesses in Figure 3.6a, and allows the GA to escape local minima. Moreover, constraining the genome length -either hard coding a maximum length or inserting a penalty term in the evaluation function where longer genomes score worse fitness valueswould end up leading the population to a local minimum. As an alternative, a periodic change in fitness goal is proposed, where during a first period of generations the algorithm seeks expression levels closer to the French flag model, and later this genome length penalty term is added with increasing weight in the evaluation function in order to obtain simpler GRNs, much in the same line as the parsimony factor used in Section 2. 3. optimization across space 35 (a) Fitness evolution across GA generations. (b) Genome length evolution across GA generations. Sudden maximum and minimum peaks correspond to gene crossover between two average length genomes yielding a very long and a very short offspring. Figure 3.6: Evolution of fitness and genome length across GA generations. (a) Obtained expression profile. (b) Comparison between normalized expression profile and theoretical expression profile (FFM: French flag model). Figure 3.7: Stable expression profiles across the spatial dimension. pr12 (green) is the protein sensing the morphogen gradient, pr00 and pr01 are the color determinants. 3. optimization across space 36 (a) Obtained 2D expression map. (b) French flag model 2D expression map. Figure 3.8: Comparison between stable expression maps across the 2D spatial domain. To obtain blue, white and red cell types, the normalized concentrations have been colormapped as follows: if the sum of normalized pr00 and normalized pr01 is smaller than 0.5, the cell type is red. If it is between 0.5and 1.5, it is blue, and if it is greater than 1.5, it is a white cell. 3. optimization across space 37 pr00 0.875 pr15 0.75 AND AND AND AND AND AND AND AND AND AND pr01 0.5 0.875 pr13 2.0 0.857 (2) 0.625 (6) AND AND pr02 0.4375 pr03 0.25 AND AND pr11 0.077 (56) 0.5625 0.667 AND pr12 0.4 (4) 1.0 AND 0.417 0.3125 AND 0.8 1.0 (7) 0.375 1.714 (5) 1.714 (5) 0.312 (7) 0.75 0.75 0.4 (4) 0.75 (3) 0.385 (24) 0.417 1.167 0.417 (18) 0.417 0.417 (6) 0.417 1.167 0.417 0.417 Figure 3.9: Final Gene Regulatory Network obtained to match response with the French flag model. Red edges indicate inhibition regulation, green edges indicate activation regulation, and black edges are inputs to AND gates or nodes self-regulation (decay rate). Decimal numbers are the apparent dissociation constant Kij (3.5), and integers between parenthesis are the Hill coefficient nof the given regulation. If not shown, n= 1. Graph layout follows an edited DOT graph description language specification [17].