Repositorio Institucional de Documentos

2014 79 José Ignacio Requeno Jarabo Formal methods applied to the analysis of phylogenies: Phylogenetic model checking Departamento Director/es Informática e Ingeniería de Sistemas Colom Piazuelo, José Manuel Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es José Ignacio Requeno Jarabo FORMAL METHODS APPLIED TO THE ANALYSIS OF PHYLOGENIES: PHYLOGENETIC MODEL CHECKING Director/es Informática e Ingeniería de Sistemas Colom Piazuelo, José Manuel Tesis Doctoral Autor 2014 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Tesis Doctoral Formal methods applied to the analysis of phylogenies: Phylogenetic Model Checking Autor: Jos´e Ignacio Requeno Director: Jos´e Manuel Colom Departamento de Inform´atica de Ingenier´ıa de Sistemas EINA, Universidad de Zaragoza We are all ignorant, but not all ignore the same things —Albert Einstein Agradecimientos Me gustar´ıa recordar que esta traves´ıa no habr´ıa llegado a buen puerto sin la colaboraci´on de diversos factores. En primer lugar, tienen un hueco importante en mi memoria las dos becas de investigaci´on (I3A [171-03] y DGA [B117/10]) y proyectos de investigaci´on del Ministerio ([TIN2008-06582-C03-02] y [TIN201127479-C04-01]) que han financiado todos estos a˜nos de esfuerzo. En segundo lugar, pero no por ello menos importante, las relaciones socio-personales tambi´en han influido de manera especial. Como es de rigor en estos casos, hago una menci´on particular a la familia, al director de tesis, a los amigos y compa˜neros del laboratorio de becarios que han cre´ıdo en m´ı y en esta tesis desde el principio, y me han apoyado en los momentos de dificultad. Como esta enumeraci´on personal suele degenerar en una lista interminable de nombres, fuente de celos, peleas y ri˜nas cari˜nosas sobre el orden de aparici´on y qui´en ha colaborado m´as en desarrollo de los acontecimientos, he optado por cortar por lo sano y evitar envidias. Por ello, a <nombre>[A-Za-z]+<\nombre>en particular, y a todo el resto de h´eroes an´onimos en general, os agradezco que hayais estado ah´ı ayud´andome y gui´andome en los momentos de debilidad. ¡Gracias! iii CONTENTS xi 10.3. Extending the Labels of the Phylogenetic Tree and the Results of ModelChecking............................. 94 10.4. Quantitative Properties in Phylogenetics . . . . . . . . . . . . . . . 97 10.5.Conclusions ............................... 98 11 Timed Transition Systems and Logics 101 11.1.Introduction...............................101 11.2. Timed Logic and Structure . . . . . . . . . . . . . . . . . . . . . . . 104 11.3. Algorithm for Timed Model Checking . . . . . . . . . . . . . . . . . 108 11.4.Conclusions ...............................111 12 Computing Distances Between Symbolic Objects 113 12.1.Introduction...............................113 12.2. Returning Time Distances as Output of the Model Checking Procedure ..................................114 12.3.Conclusions ...............................118 13 Conclusions 119 V Quantitative Extensions of Kripke Structures and Logics for the Analysis of Phylogenies: Approaching to Probabilistic Properties 121 14 Discrete Time Probabilistic Transition Systems and Logics 125 14.1.Introduction...............................125 14.2. Discrete time Probabilistic Logic and Structure . . . . . . . . . . . 126 14.3. Algorithm for PCTL Model Checking . . . . . . . . . . . . . . . . . 129 14.4. Model Checking Tools and Experimentation . . . . . . . . . . . . . 130 14.5.Conclusions ...............................132 15 Continuous Time Probabilistic Transition Systems and Logics 137 15.1.Introduction...............................137 15.2. Models of DNA Evolution . . . . . . . . . . . . . . . . . . . . . . . 138 15.3. Maximum Likelihood Estimation . . . . . . . . . . . . . . . . . . . 140 15.4. Continuous time Probabilistic Logic and Structure . . . . . . . . . . 142 15.5. Algorithm for CSLTA Model Checking . . . . . . . . . . . . . . . . . 146 15.6. Model Checking Tools and Experimentation . . . . . . . . . . . . . 148 15.7.Conclusions ...............................155 16 Conclusions 157 xii CONTENTS VI Conclusions and Further Remarks 159 17 Conclusions 161 18 Future Work: Parametric Temporal Logic 165 18.1.Introduction...............................165 18.2. Parameters in boolean model checking . . . . . . . . . . . . . . . . 166 18.3. Parameters in timed model checking . . . . . . . . . . . . . . . . . 167 18.4. Parameters in probabilistic model checking . . . . . . . . . . . . . . 167 Bibliography 169 List of Figures 2.1. Translation from a phylogenetic tree to a Kripke structure. . . . . . . . 19 2.2. Evaluation of temporal logic operators. . . . . . . . . . . . . . . . . . . 21 4.1. Boxed nodes indicate the states where the corresponding property is held...................................... 40 5.1. Boxed nodes indicate the states where the corresponding property is held...................................... 50 5.2. Mapping of the phylogenetic tree of Figure 2.1 in SMV. . . . . . . . . . 55 5.3. (a) Time is linear with respect to set size and (b) quadratic with respect tosequencelength. ............................ 56 5.4. Counterexample of a back mutation property. . . . . . . . . . . . . . . 57 7.1. Distributed verification of EXφthrough the parallel execution of φin thedirectsubtrees.............................. 70 7.2. Time required in NuSMV for the initialization of a phylogenetic tree with GenBank identifiers in the nodes. . . . . . . . . . . . . . . . . . . 71 7.3. Memory required in NuSMV for the initialization of a phylogenetic tree with GenBank identifiers in the nodes. . . . . . . . . . . . . . . . . . . 72 7.4. Division of the original Kripke structure into two slices. . . . . . . . . . 74 8.1. Workflow diagram with the alignment, phylogenetic tree and properties asinput.................................... 82 10.1. Phylogenetic tree for the Hominoidea. . . . . . . . . . . . . . . . . . . 91 10.2. Phylogenetic tree labeled with quantitative information. . . . . . . . . 96 11.1. Phylogenetic tree and its transition system labeled with time intervals in the branches and taxon identifiers in the nodes (the DNA sequences are omitted for readability). . . . . . . . . . . . . . . . . . . . . . . . . 108 12.1. Description of a phylogenetic tree in PRISM syntax. . . . . . . . . . . 117 xiii xiv List of Figures 14.1. Mapping of the phylogenetic tree of Figure 2.1 in PRISM. . . . . . . . 134 14.2. Time required for the verification of a set of probabilistic formulas with respect to the number of tips in the phylogeny. . . . . . . . . . . . . . 135 15.1. Model of DNA substitution. . . . . . . . . . . . . . . . . . . . . . . . . 139 15.2. Unfolding of a model of DNA substitution. . . . . . . . . . . . . . . . . 140 15.3. Description of the Jukes-Cantor model in PRISM syntax. . . . . . . . . 149 15.4. Representation in PRISM syntax of the MLE equations for Figure 2.1. 150 15.5. Time required in PRISM for the evaluation of the maximum likelihood equations in a binary phylogenetic tree. . . . . . . . . . . . . . . . . . . 152 15.6. Time required in PRISMopt for the evaluation of the maximum likelihood equations in a binary phylogenetic tree. . . . . . . . . . . . . . . 153 15.7. Rewriting of the MLE equations for Figure 2.1. . . . . . . . . . . . . . 154 15.8. Time required in PRISM for the computation of the upper bound of the maximum likelihood value in a binary phylogenetic tree. . . . . . . 155 List of Tables 4.1. (1) Summary of the most important phylogenetic properties (Type T: tree; S: sequence; Q: quantitative). . . . . . . . . . . . . . . . . . . . . 41 4.2. (2) Summary of the most important phylogenetic properties (Type T: tree; S: sequence; Q: quantitative). . . . . . . . . . . . . . . . . . . . . 42 5.1. List of some available model checkers. . . . . . . . . . . . . . . . . . . . 49 5.2. Seconds needed for the creation of the Kripke structure and the storage ofproteinsequences. ............................ 52 5.3. Megabytes needed for the creation of the Kripke structure and the storage of protein sequences. . . . . . . . . . . . . . . . . . . . . . . . . 52 10.1. Summary of the most important phylogenetic properties (Type B: boolean; N: numeric; P: parametric). . . . . . . . . . . . . . . . . . . . 99 xv Part I Introduction and Formal Definition of the Framework Chapter 1 Introduction “You are to bring into the ark two of all living creatures, male and female, to keep them alive with you. Two of every kind of bird, of every kind of animal and of every kind of creature that moves along the ground will come to you to be kept alive. You are to take every kind of food that is to be eaten and store it away as food for you and for them.” Noah did everything just as God commanded him. Genesis 6:19-22 1.1. Biological Context: Phylogenies and Phylogenetic Analysis The study and exploration of life is, for reasons as evident as they are difficult to formulate, one of the most deeply rooted constants in human thought. Understanding and classification are intimately related, and consequently it should come as no surprise that the branches of science that concern themselves with biological classification figure prominently in the history of modern science, from the revolutionary breakthroughs of Linnaeus and Darwin to the newer disciplines that have evolved in their wake. Generally, the goal of taxonomy, into which Linnaeus breathed new life in the 18th century [113, 114, 115, 116], is to define and organize populations of organisms as taxonomic units (taxa) in hierarchies according to shared traits; in the domain of life, these may correlate in varying degrees to relations of descent. Additionally, systematics is informed by Darwin’s 19th-century defense of common ancestor [54], by now universally accepted by the scientific community, and therefore seeks to propose phylogenies that reflect said existent descent relationships. Regardless 3 4CHAPTER 1. INTRODUCTION of their finer points, an overwhelming majority of taxonomic and systematics approaches commonly agree in their basic conception of evolution and classifications as hierarchical trees, whether their splits represent speciation events (mutations or sexual reproduction) and groups are founded on relevant attributes. Phylogenetic trees, or phylogenies for short, continue to be useful abstractions for modeling and characterizing evolution over time [4]. Biologists consider them as feasible, traditional and wide-spread approaches of more powerful data structures such as phylogenetic networks or pedigrees. It is accepted that neither inheritance is single and infallible, nor descent strictly linear, with sexual reproduction and horizontal transfers of genetic material blurring this ideal model to some extent. Yet the lines are clearly visible, and the tree model is perfectly valid for vast regions of the postulated tree of life, if only we shift our focus from individuals to populations and beyond [18]. Biologists use phylogenetic trees in many different ways to solve both scientific and practical problems. For example, they use trees to make predictions about fossils or poorly-studied species, and to learn about the evolution of complex features, the order of evolution or the evolution of diversity [167]. During the last century, systematics and taxonomy have provided a wide range of methodologies, both theoretical and technical, which have contributed to a better understanding of evolution and the building of the tree of life. Among the former, Hennig’s 20th-century cladistics [89] has proven a solid and durable theoretical methodology, aided in recent times by the development of computational phylogenetics. In addition, the concept of evolutionary distance allows for exploring tangible numerical relationships between sets of populations or individuals characterized by biological features such as DNA. Computer science tools have upgraded the capabilities of biologists for the tree construction. A great amount of software packages for tree-building guides the power of computer science towards distance-based [148] (PHYLYP [70]) and character-based methods [74] (PAUP [161]) or maximum likelihood estimations [69, 71, 176] (RaXML [155], jModelTest [53]). In particular, the maximum likelihood estimations contribute with statistically consistent techniques for inferring the probability that the phylogeny accommodates a specific evolution drift, provided by the scientists in terms of DNA substitution models [117]. On these days, one of the most relevant challenges orbits around the validation of the inferred tree. Both disciplines, taxonomy and systematics, make heavy use of empirical data, with proposition, verification and generalization of hypotheses over the reconstructed tree playing a central role in their application [130]. Once trees are inferred, the objective is to extract or to discover the implicit biological messages and structural properties underlying in the phylogeny under inspection. For in- 1.4. OBJECTIVES OF THIS PHD THESIS 11 or conserved regions in a clade. Moreover, the idea of hypothesis verification works in a closed-loop, helping the phylogeneticist to refine or discover properties using counterexamples obtained from unfulfilled properties. To this end, the introduction of model checking techniques is summarized under three different perspectives: 1. First of all, our goal is to separate the implementation (software tools) from the analysis of a phylogeny, denoted with a set of properties the user desires to check. The objective is to develop a different kind of phylogenetic analysis. Inside a model checking framework, the analysis and study of properties is realized independently of the software tool that makes the particular computations. That is, the model checker tool is a generic package that interprets the specification of a system and automatically realizes the computations to carry out the analysis. Then, the biologists can concentrate their activity in the specification and interpretation of the results. In addition, the model and specifications are usually defined as separated input files for the model checking tool, allowing the verification of the same properties over multiple trees and the interpretation of different results or behaviors. This is a very instrumental reason, and it sometimes presents a conflict between generality and efficiency. Nevertheless, there are many developments inside the formal verification domain that are transparent to the application domain and that can be translated to the context of phylogeny. Therefore, there is no need to return to the path that the computer scientists already visited before. For example, the research and implementation of methods for the distribution of computations in order to speed up the performance, or the extraction of efficient representations and symbolic manipulation of the information. The use of experienced software with up-to-date technologies and a good community support such as NuSMV [45] facilitates the study of phylogenetic properties. Our objective is to take advantage of these tools and adapt them to the domain of phylogenetics if necessary, instead of coding a completely new software or focusing on particular algorithms. On the other hand, current packages for phylogenetic analysis such as Mesquite [120], ape [67] or Pycogent [105] are becoming really complex with thousands of code lines. They are difficult to maintain, even when they include specific and adapted algorithms for this field. The attractiveness of model checking comes from the fact that it is completely automatic, generic and independent of the model and specifications (i.e., the learning curve for a user is very gentle). 2. The second cornerstone of our approach is the introduction of the symbolic manipulation of biological objects used to reason inside the world of phylogeny. In this sense, the introduction of formal definitions through logical 12 CHAPTER 1. INTRODUCTION formulas, interpreted over the phylogenetic tree, allows to characterize sets that are not explicitly enumerated. These sets of objects are symbolically managed and it is possible to define specific calculus to manipulate them. In a complex world with many levels of abstraction, such as the phylogeny, the ability to systematically manipulate components or objects represented in a symbolic way constitutes a key advantage. This requires some abstraction tool allowing reasoning with the right objects at each level, removing unnecessary details for the searched results. For example if we are trying to work in the field of cladistics, then one of the right objects to reason with is the set of clades. Hence, we need to have a symbolic representation of these objects and a calculus to reason with them. 3. Last but not least, we propose model checking for studying complex systems with single or multiple phylogenetic trees as input. Our objective is to obtain and to organize the information contained in these structural-logical models. Here we try to extract properties of a phylogeny, combine them, compare the results among trees and split the properties into smaller meaningful ones. In other words, we need an algebra to manipulate properties symbolically. This is a difference with respect to the classical verification framework, which restricts the analysis to the verification if a property is satisfied or not. We propose a classification of phylogenetic properties and a methodology for the specification of these hypothesis using temporal logics. In addition, our phylogenetic model is initially represented by a phylogenetic tree labeled with the genome of each taxon, but the states of the phylogeny can be also enriched with extra information complementary to that obtained from the DNA sequences or the output results obtained from the evaluation of previous properties, for example, with phenotypic information. Therefore, the model checking techniques work with properties that can be interpreted over different models (phylogenetic trees) or over phylogenies that maintains the structure but the information content is changed. We also work with models that include probabilities (Markov chains). Modifying the behavior of the logics and data structures, we can work in the future with phylogenetic networks or pedigrees as well. To sum up, this approach allows the specification of properties in temporal logic in an unambiguous and mathematically strict way. Hence, biologists can formalize consensus definitions for controversial phylogenetic terms, mainly in quantitative analysis. The definition and analysis of a new phylogenetic property doesn’t need the implementation of any new algorithm, which reduces the risk of coding errors and simplifies the software cycle of life. The same property can be exported to and tested in different models (phylogenetic trees) with minor changes, fixing them 1.5. ORGANIZATION AND LIST OF PUBLICATIONS 13 for phylogenetic networks if necessary. It also offers counterexamples in case of a model fails to satisfy a property serving as indispensable for the debugging. The previous arguments correspond to the conceptual advantages of model checking beyond the computational efficiency. In fact, we try to demonstrate that our contribution consists of the adaptation of the model checking framework as an unifying formalism to realize many different tasks under the same generic formal description languages and tools, offering the possibility to create different abstraction levels throughout the symbolic creation and management of objects. The work introduced in the first theoretical sections is completed with a second part oriented to the usability and competitiveness of our approach. We have evaluated the performance of existing model checking tools in the field of phylogenetics. Due to the particularities of the biological context and the huge volume of data they have to manipulate, we have discovered that most of the software packages suffer from slowness. Then, we are forced to design, adapt, propose and implement new solutions for increasing the efficiency of our methodology in phylogenetics. Finally, the formalism is enriched with time and probabilities for model checking. 1.5. Organization and List of Publications The contents presented in this thesis is supported by several publications in international conferences and journals. The introduction of model checking as a formal framework for analyzing biological properties led us to our first seminal paper [26]. In [138], we develop the verification of a particular phylogenetic example with proteins coded by human mitochondrial DNA (mtDNA). There, we used a generic and automated model checking tool. The journal paper [143] summarizes the basics of previous papers and extends the operational of temporal logics to new phylogenetic properties unconsidered until now. Due to the peculiarities of the phylogenetic analysis (i.e., mainly huge trees with a considerable amount of information in the states, such as the DNA), the model checking tools need to be optimized in order to keep a reasonable performance and ease the state explosion problem. The paper [139] presents the advantages of slicing the DNA strings associated to the states, and the partition of the phylogeny into small subtrees. Besides, the original idea of sliced model checking is complemented with the addition of an external database for storing the sequences. The journal paper [140] gathers the notions introduced in [139] together with the implementation and preliminary results of [141]. At the end, the thesis extends the standard temporal logics with explicit time and probabilities so as to manipulate and interrogate the tree with quantitative 14 CHAPTER 1. INTRODUCTION information. The paper [142] exemplifies the necessity of introducing probabilities and discrete time for the phylogenetic analysis of a real phenotype. Therefore, the thesis is arranged in six parts. Following this introduction, Part I summarizes the essentials of phylogenetics and model checking. It explains the roots which bridge model checking and phylogenetic analysis: phylogenies as logical models, phylogenetic specifications as temporal logic formulas, and automated system verification via model checking tools. Later, Part II describes the logical specification of non-trivial structural properties of cladistic classification and sequence composition. It also details the key steps for implementing phylogenetic trees and biological properties within the scope of a Symbolic Model Verifier (SMV) tool in order to obtain both feasibility and performance criteria for our approach. Part III considers different compatible optimizations for scaling the framework for bigger systems, mainly the introduction of tree partitions, state slicing and external databases. A workflow describes the interconnection and cooperation of the modules implementing the mentioned techniques. Next, the logic and structures presented in Part II are extended for quantitative purposes (explicit time and probabilities) in Part IV–V. With a set of real examples, we motivate the interest of inserting explicit clocks and likelihoods in phylogenies. Finally, Part VI gathers the conclusions drawn from this research and outlines the future work. We highlight that the introduction of parameters in models and specifications allows the symbolic manipulation of relations, converging to a potential automatic discovery of properties and model exploration. Chapter 2 Bridging Worlds: Phylogenetics and Model Checking The fact is, what we’re doing could be construed as, forgive me sir, collaboration with the enemy. Perhaps even as treasonable activity. . . Must we work so well. Must we build them a better bridge than they could have built for themselves? —Maj. Clipton, The Bridge on the River Kwai 2.1. Introduction Transition systems are considered powerful formal models for the study of concurrent systems [13, Definition 2.1]. They are machines composed of a set of (probably) infinite states labeled with information describing the system behavior at that point, and a set transitions among them. Formally, a state transition system is a pair (S, R) where Sis a set of states and R⊆S×Sis a binary relation (transitions) over S. Given a model of a system and a suitable logic to reason about its behavior over time, it is possible to achieve automated, exhaustive verification of properties of interest. In this section we show that phylogenies can be understood and represented by such models. Phylogenetics and model checking can be bridged after reflecting on some considerations about the processes of modeling and specification. We begin by identifying the foundations of evolution with those of transition systems and continue 15 16 CHAPTER 2. BRIDGING WORLDS with the study of their temporal nature and their logical formulation. We conclude with an overview of verification under the model checking paradigm. 2.2. Evolution as a Transition System At the highest level of abstraction, phylogenies postulate partial models of the evolution history of sets of living organisms. They can be represented as directed graphs, though phylogenetic trees are widespread representations and they suffice for most common purposes; phylogenetic trees will be adopted in this thesis. Rooted labeled trees offer a realistic model of aggregated evolution, in which each vertex represents an inferred state of the evolution characterized by a population of related individuals who mate among themselves and are denoted by a common, compatible heritage (e.g., biological sequences such as DNA) or other information such as their attributes (e.g., morphological and physiological data). The vertex are arranged along the paths of the tree according to a evolutive process. The methodology presented here is extensible to other domains (trees or networks) and phylogenies annotated with different kind of information. Although phylogenetic trees don’t match exactly with a gene tree [131, 132, 119], we start focusing on trees built from genes and genomic alignments for simplicity. Transformative events that modify heritable information give rise to new states which are reflected in oriented parent-child edges in the graph. Transitions are comparable to “instantaneous” speciation events or mutations. Note that neither explicit time nor the ordering of child states are part of the model, which is consistent with the descent semantics of the tree. In this context, time is implicitly represented as a relation order that indicates the direction of the evolution drift caused by the speciation events. Definition 1 (Rooted Labeled Tree).Let Σbe a finite alphabet and la natural number. A phylogenetic tree over Σlis a tuple P= (T, r, D), where: T= (V, E)is a tree graph, r∈Vis its root, and D:V→Σlis a dictionary function that labels each vertex with its associated taxon information, mainly the genome sequence. Gene trees are typically built from finite sets of words of uniform length resulting from alignment algorithms. These words, which commonly represent presenttime taxa, are required to be in bijective correspondence with the set of leaves of the trees. The inclusion of extinct species, which are improbably direct ancestors of extant taxa, are also located in terminal nodes. Ancient DNA sequences far 2.2. EVOLUTION AS A TRANSITION SYSTEM 17 before 100,000 years are unusual in the phylogeny due to its deterioration, but they can be inferred for the internal nodes of the tree using maximum parsimony for the direct common ancestors. Restrictions can be added or removed to adjust the phylogeny as needed (i.e., including horizontal transference, sexual reproduction or multiple roots for phylogenetic networks), but in any case the nature of trees as reactive systems should become clear by now. They are composed of independent states of evolution that interact indefinitely with their environment. One of the most prominent features is their state, i.e., their hereditary information or a suitable portion thereof. Consequently, it is possible to naturally effect modeling and verification of evolutionary systems, and to this end data structures for the representation of transition systems become a very attractive solution. A Kripke structure provides a graph-based data arrangement that represents the behavior of the system and provides semantic information for querying properties in temporal logics [122]. Definition 2 (Kripke Structure).Let AP be a set of atomic propositions, i.e., boolean predicates that describe the observable properties of a state. A Kripke structure over AP is a finite transition system represented by a tuple M= (S, S0, R, L, AP), where: Sis a finite set of states, S0⊆Sis the set of initial states, R⊆S×Sis a total transition relation between states, i.e., for every state si∈S, there exists sj∈Ssuch that (si, sj)∈R, and L:S→2AP is the labeling function that associates each state with the subset of atomic propositions that are true of it. A Kripke structure models a system that is capable of an infinite number of behaviors or paths, infinite sequences of successive states π=s0s1s2. . . such that s0∈S0and (si, si+1)∈R, i ∈N. The set of possible executions (paths) in a structure can be unfolded into its computation tree1. Here we focus on Kripke structures that represent a certain tree as a (hopefully real) computation of the evolutionary process, or rather the set of computations that result in the hypothesized patterns of evolution. Relations between states and atomic propositions, and between tree branches and state transitions, are of utmost in importance, and raise some interesting issues that need to be addressed: 1A computation tree is a rooted tree-like structure in which the future is not determined by the current node and, therefore, there are many different paths reachable from this state. 18 CHAPTER 2. BRIDGING WORLDS Ideally, the state of an evolution process is uniquely identified by its biological information. Sequences also determine the atomic propositions that form the basis of logical properties: the presence of a certain alphabet symbol at a given position. If separate vertices of the tree must share a sequence, their states can be enriched with auxiliary properties to preserve the unique identity of each. Once a one-to-one correspondence between tree vertices and states has been established, it remains to resolve the identification of branches as state transitions. Transitions are comparable to “instantaneous” speciation events. In sum, trees are essentially a present-time, local, tentative snapshot in the execution of a potentially infinite evolution system. In order to translate trees to the infinite-path semantics of Kripke structures, we need to add self-loops to terminal vertices so as to deadlock their states. At this point, we can define a suitable branching-time structure for phylogenetic trees for the interpretation of temporal logic formulas which express properties in them. The most common state formula determines whether the present state is associated with a sequence seq =σ1σ2. . . σl∈Σl(or possibly a partial sequence or a set of sequences). Sequences will be manipulated symbolically, as will sets, as the aggregation of their parts. In logical terms: seq ≡Vl i=1seq [i] = σi(2.1) These expressions related with the DNA compose the basic set of atomic propositions (AP) that are tagged to the states, that is, AP ={seq [i] = σ|σ∈Σ, i≤l}. The set of AP accepts the inclusion of new information apart from the genome. Definition 3 (Branching-time Phylogeny).A tree (per Definition 1) P= (T, r, D) is univocally defined by the Kripke structure M= (V, {r}, R, L), where: Ris the transition relation composed of the set of tree edges (directed from r) plus self-loops on the leaves: R=E∪ {(v, v) : @(v, w)∈E∧v, w ∈V} with v6=w. Lis the standard labeling function defined by AP, under which a state v mapped to D(v) = seq with seq =σ1σ2. . . σlsatisfies the family of properties seq [i] = σi,1≤i≤l, plus a unique state identifier in the case of several states sharing the same atomic propositions. Roughly speaking, we can say that the Kripke structure reflects the parentchild relations in the original phylogenetic tree representing the evolution process. Self-loops in terminal nodes allow for infinite computation paths in the original 2.3. TEMPORAL LOGIC AS A SPECIFICATION LANGUAGE 19 Figure 2.1: Translation from a phylogenetic tree to a Kripke structure. phylogenetic tree. The labeling of states with atomic propositions in the transition system helps in the identification of nodes during the verification process. Thanks to the codification explained, the detection of a particular instance of a DNA sequence is translated to the selection of states satisfying a set of logical equalities with the form seq [i] = σi,1≤i≤l. In this way, we can manipulate sequences symbolically. In case of collision, a explicit definition of state identifiers ensures uniqueness. The set of atomic propositions has to be able to represent distinct populations that share a common sequence distinguished as different states. Figure 2.1 illustrates this translation process from a (phylogenetic) rooted labeled tree (Definition 1) to a Kripke structure (Definition 2). It should be emphasized that our proposal can be used in future works for other phylogenetic structures such as phylogenetic networks [52] by reformulating Definitions 1–3 as required. Although they are more generic and powerful representations of the evolution process, the complexity for interpreting and manipulating specifications in temporal logic in a phylogenetic network places this data structure out of our scope in this thesis. 2.3. Temporal Logic as a Specification Language Temporal logics are formal systems that allow the representation and manipulation of logical propositions qualified in terms of time [134, 122, 46, 63, 13]. In 20 CHAPTER 2. BRIDGING WORLDS the context of transition systems, they are used to define properties on sequences of transitions between states of a system through a convenient abstraction of it (in the present case, a specific kind of Kripke structure). For example, properties may express whether it is possible that one particular type of state may be reached at a particular point or whether a certain property will always hold. Temporal logics can be classified according to how they treat sequences of events: whereas linear-time logics (LTL) [134] deal with individual paths π= s0s1s2. . .,branching-time logics take into account the set of possible progressions from each state, hence reasoning globally about the computation tree. Computational Tree Logic (CTL) [46], a versatile exponent of the latter, has been widely adopted by the model checking community. Phylogenies represent evolutionary processes that are mainly branching in nature due to the (hypothetically independent) speciation events over the time. Branching-time logics in general, and CTL in particular, are remarkably well suited to the description of phylogenies because the structure of phylogenetic trees matches with the topology of a computation tree and CTL provides a set of operators for inspecting and manipulating those paths. We propose Phylogenetic Tree Logic (PTL) as a temporal logic close to CTL for the specification of evolutionary processes. PTL reinterprets the quantifiers of first-order logic as path quantifiers, expressing the fulfillment of a property throughout all computation paths (A), or at least one computation path (E). These two must be immediately qualified by one of five temporal operators, of which three express the satisfaction of a property eventually in time (F), at all times (G), or in the next state (X); and two are conditional constructs in which a precedent is verified until a consequent comes into force (U), or until and including the moment when it does, if it does (R). A complete grammar and semantics of PTL formulas can be defined from a minimal and complete subset of logical operators (see Figure 2.2). Definition 4 (Phylogenetic Tree Logic).A temporal logic formula φis defined by the following grammar, where p∈AP: φ::= true |p| ¬φ|φ∨φ|EX (φ)|EG (φ)|E[φUφ] (2.2) The formulas are checked against a structure Mconsidering all paths πfrom a certain state s0. Notice that M, s0φmeans that s0satisfies φ. The semantics of well-formed formulas is as follows (let π=s0s1s2. . .): M, s0p⇔p∈L(s0), M, s0¬φ⇔M, s02φ, M, s0φ∨ψ⇔M, s0φor M, s0ψ, 3.1. FOUNDATIONS OF THIS THESIS 27 Probably, they alert with counterexamples that report erroneous relations in the phylogeny. The result of testing phylogenetic properties expressed in temporal logics may confirm or refute biological hypothesis about the speciation drift, which assesses the tree quality with respect to different criteria. For example, rootedness is an important factor in practice and should not be taken for granted, even for simple phylogenetic trees. An overwhelming majority of phylogeny reconstruction methods, including all those in common use, produce as their results undirected trees, which lack a distinct root unless ascertained by external means, e.g. by using outroots. To this end, model checking is useful as a method to locate the root and discriminate among several potential candidates. Given a set of properties, model checking will consist of detecting the set of potential roots that verify them. That is, consider a property to verify on an unrooted phylogeny. For instance, we want to inspect the roots whose mean distance to the leaves is over a certain threshold, i.e., is the tree balanced if we take that node as root? (see Chapter 12.2). It is possible to determine the set of nodes which can be the root of the tree according to that property by attempting to verify the property taking each of the nodes of the tree as its root. This involves the symbolic manipulation of sets of states by means of the propositions they satisfy. The total cost of verification is linear with respect to the number of potential roots, which nonetheless is composed of completely independent problems (though they may have subproblems in common, etc.). This can be parallelized with perfect speed up. Each property in a collection restricts the possible locations of the root to a subset of the nodes of the tree. We may conclude that a node is consistent with the rootedness of the tree under all those properties if it is acceptable as a root for each and every one of the properties. The number and type of properties validated over the tree gives a way to judge the strength of the root under consideration. The feedback returned by the model checking process helps to update the phylogenetic structure. In spite of the mentioned advantages for defining and evaluating phylogenetic hypothesis, the model checking techniques has never been introduced in the domain of phylogenetics yet. They are mainly employed for industrial domains and particular operations. Model checking is barely introduced in computational biology, where it has been applied to fields in which (mostly) quantitative properties over temporal data are analyzed. This fact contrasts with the principal limitations of phylogenetic analysis. The goals of introducing model checking techniques in phylogenetics have been a) the presentation of a formal framework for the investigation of biological hypothesis that solves the lack of independence between specification and verification, and b) the use of generic and powerful tools for the validation of those properties and the symbolic manipulation of states instead of 28 CHAPTER 3. CONCLUSIONS the codification of new software. As final remarks, temporal logics assume a well-defined flow of time, and thus are dependent on an oriented graph. Kripke structures are oriented by definition, i.e., there exists a initial state at least. Branching-time temporal logics define properties over the computation tree of a Kripke structure. A rooted phylogenetic tree matches with the notion of a computation tree because the speciation events are similar to the transitions and the taxa correspond to the states of the computation tree. In case of unrooted phylogenies, it is necessary to either adapt a temporal logic to accept uncertainty in the flow of time, or the model checking algorithm. Clearly, the Kripke structure needs to be generalized. In any case, it is possible to overcome this problem by defining a partial orientation in regions where the evolution flow is obviously known. In sum, the phylogenetic tree is still a feasible and useful model over which we inspect biological hypothesis despite the limitations and drawbacks it has. The inherent hierarchy of the data structure allows the adaptation and optimization of the model checking techniques to this domain. Besides, the use of a phylogeny as a model of the evolution contrasts with the use of DNA substitution models, which are considered as the default models for analyzing the principles operating behind the macro evolution process. More complex and refined data structures like phylogenetic netwoks can be considered in the future. 3.2. Hypothesis: Phylogenetic Trees versus Phylogenetic Networks As well as phylogenetic trees, other structures such as phylogenetic networks are suitable for phylogenetic analysis because they represent a generalization over trees [52]. In particular, they exhibit an extra level of abstraction for modeling evolution based on the inclusion of potential multiple roots, horizontal transferences, sexual reproduction and the existence of cycles. In fact, the inclusion of cycles in the phylogenetic network is an useful idea. This involves the addition of new semantic information to the phylogeny. For example, a phylogenetic network can be understood as a compact representation of a set of phylogenetic trees where branches and nodes overlap in the same structure [172]. Thus, the verification process focuses not only on the validity of the logical proposition but the detection of trees in which they hold. The identification of which nodes belong to each tree is easily solved by defining the membership with a explicit variable in every node. This feature represents an efficient way of processing several trees together. Cycles in a phylogenetic network could also stand for the horizontal transference 3.2. PHYLOGENETIC TREES VERSUS PHYLOGENETIC NETWORKS 29 of fragments of sequences (genes) between individuals of simple organisms like bacteria, while they would represent sexual reproduction for the rest of animals and plants. In this case, there would be only one phylogeny under study but the complexity of the verification process would be increased due to the explosion of paths. As the notion of cycles is already included in the definition of Kripke structures, and both phylogenetic interpretations of cycles are compatible with the semantics of temporal logic, phylogenetic networks are also feasible in the context of model checking. Relevant phylogenetic properties expressed in temporal logics can be imported to phylogenetic networks with fairly minor changes. Then, phylogenetic networks are more powerful and expressive. However, we maintain the use of phylogenetic trees in the rest of this dissertation for simplicity and compatibility with previous studies: phylogenetic trees continue to be useful and widespread abstractions of the evolution over time. This thesis places the essentials for the introduction of the verification process of phylogenetic properties using model checking, and therefore trees suffices for the first steps. Besides, a tree-like structure simplifies the use and adaptation of model checking. In any case, we desire to remark that the work developed here is flexible enough for the extension of phylogenetic analysis to other data structures in the future, such as networks. Part II Phylogenetic Properties: A Logic for Qualitative Logical Properties The inspection of properties over a phylogeny needs the adoption of a temporal logic capable of capturing the essence of these biological hypothesis. The specification of phylogenetic properties with temporal logic is presented here for the first time. The existence of coalescence, conservation or covariation of sequences, gene inversions, reversions and a complete set of mutation patterns can be analyzed with a basic temporal logic operating over the paths of the tree. They allow the study of the topological structure, obtaining meaningful information about the peculiar nucleic arrangements and aggregations over the tree, and possibly examining the phylogenetic quality by recognizing anomalous characteristics. The first chapter of this part is dedicated to classify the kind of properties that can be represented with a branching-time temporal logic and illustrate the methodology for its definition. By default, we organize the properties in three families depending on the phylogenetic information they require: a) topological or tree properties, that need the structural information of the phylogeny, b) sequence or state properties, that work with the sequences attached to the nodes, or c) a combination of both. For each case, we proceed incrementally with a modular composition and specification of subproperties in order to show every step of the design. Subsequently, the next step faces the verification of those phylogenetic properties over a specific model checking tool. There, we show how to translate the phylogenetic tree into the syntax of a generic software and code the specifications accordingly. In case of failure to accomplish the property, a counterexample is returned by the tool. We use an example to clarify the interpretation of the counterexample output. Finally, the performance results demonstrate the feasibility of our approach over human mitochondrial DNA data and sketches future work. 33 Chapter 4 A Temporal Logic for Phylogenetic Analysis No, no, you’re not thinking; you’re just being logical. —Niels Bohr 4.1. Introduction This chapter describes a methodology to specify biological properties of phylogenies employing the previous logical framework. Syntactically, a non-trivial property can be broken down into simpler, meaningful ones, and synthesized from these. The benefits of such logical decompositions are twofold: first, they simplify formalization and favor readability; second, modular properties can be reused in complex constructs, and variations on a formula can be produced by local adjustments. The main idea is to apply formulas in temporal logic to extend the aptitudes of current phylogenetic analysis tools and cover the increasing necessities of biologist for the verification of hypothesis (see Section 1.3). These tasks include the detection of common ancestors, monophyletic clades [64], reversions (back mutation) [65, 1], and probably extend the hypothesis testing with quantitative information for the study of speciation rates (extinction and diversification) and mutation rates in DNA substitution models [117, 53]. Four classes of queries can be identified: global,topological or tree properties, in which the structure of the phylogeny itself is placed under scrutiny and relations between taxa inspected; local or sequence properties, where compositional sequence features take center stage, possibly aided by the placement of constraints; a combination of both, where the tree topology and the sequence alignment are 35 36 CHAPTER 4. TEMPORAL LOGIC simultaneously required; and quantitative properties, that extends the preceding properties with explicit time and probabilities. In the following sections we describe the relevant properties of each group and exemplify their formal modeling (quantitative properties are studied apart in Part IV–V). Table 4.1 summarizes relevant phylogenetic properties according to the classification introduced. 4.2. Tree Properties Tree properties are mainly of a cladistic nature, asking about the organization of the phylogenetic tree in compatible regions and the biological characteristics that are satisfied in them [89]. One of the most frequent basic queries is whether a set of extant organisms Sunder study constitutes a monophyletic group or clade in a phylogeny. That is, does the phylogenetic tree contain a subtree that has exactly those organisms as its leaves? Formally, the PTL formula over the Kripke structure will be true if there exists somewhere in the tree a reachable state (EF) which is the root of the group S, and thus: a) everything that has to be in, is in; and b) everything that has to be out, is out (there are no outsiders, or conversely, only insiders can be found inside): clade (S)≡EF (in (S)∧out (S)) (4.1) The inclusion rule states that for each unique sequence sof the set S(and-logic V) there is a path that finds it as its leaf; the exclusion rule demands that all paths end in a leaf from the same set. Within infinite computations, leaves are found either through a terminal boolean variable or through the pattern F AG (s), since by construction whole uniform computation subtrees can only be found precisely as a consequence of leaf states: in (S)≡Vs∈SEF AG (s) (4.2) out (S)≡AF AG Ws∈Ss(4.3) As noted, individual properties in isolation have useful semantics in their own right. Here, in (S) is satisfied by all containing clades, and out (S) by all strict subclades. In particular, the roots (non-terminal nodes) satisfying in (S) define the group of common ancestors (CA): CA(S)≡ ¬terminal ∧in(S) (4.4) By default, the set of taxa, and consequently their sequences, are located in the tips. However, if the phylogeny is extended with both ancestral and leaf sequences, the structure of the clade property remains unaltered. Only the inclusion and 4.3. SEQUENCE PROPERTIES 43 As a common understanding, closely related members share analogous genomes, satisfying identical patterns and properties. For example, the preservation of biochemical properties is essential for maintaining the stability of the resulting protein and the viability of an organism, placing important constrains to the evolution [86]. The amino acids conform the raw material of the protein, and consequently they determine the three-dimensional structure and singularities of the molecule. Proteins develop crucial functions in living organisms. The modification of a nucleic base, and indirectly its associated amino acid, may bring a recalibration in the molecular morphology which damages its functionality and propagates the instability to a chain of biological reactions [50, 173]. Therefore, natural selection favors the persistence of the molecular structure limiting the mutations to compatible amino acids. The existence of synonymous substitutions2provides some flexibility to the change of nucleic bases, but in any case the conservation of biochemical properties and sequences reaffirms the use of patterns. When a nucleic modification inevitably changes the protein morphology or alters the biochemical equilibrium, the effects in the organism are unpredictable. In the worst case, the mutation becomes lethal and punishes the survival of the individual. In the best situation, it is ubiquitous or a new evolutionary advantage arises at the conclusion of this process. In an intermediate context, a single mutation occasionally triggers a group of (intra or inter-gene) changes that neutralize the consequences of the original modification in a short period of time. In other words, they try to compensate the attributes of the new amino acid in order to maintain the biochemical equilibrium. This character coevolution is often captured in the vicinity of the phylogeny [77]. The number and degree of pattern violations (nucleotide and biochemical conservation, protein shape) together with their circumstances assists to the discernment of the repercussions caused by the mutation. According to the previous comments, patterns represents global correctness constraints that are supposed to hold across the whole phylogeny. They can be categorized as follows: Conservation is modeled as a restriction on the symbols that can occur at a given position in a sequence. Commonly, the pattern will codify a unidimensional boolean table that classifies each symbol as permissible or impermissible. However, it is possible to define general families of compatible elements, not bounded to specific positions, and restrict their usage to exactly one of these positions, among other extensions. Covariation imposes a relation of dependence between two (or more) positions in a sequence. It can be represented as a bidimensional boolean table 2Codons that transcript the same amino acid 44 CHAPTER 4. TEMPORAL LOGIC which states, for each symbol in the first column, the set of symbols that may appear in the second column. Typically, for the property to be meaningful, associations between symbols will be sparse. A combination of both. A global pattern thus defined is easily verified by extending it over the computation tree: global (p)≡AG (p) (4.16) Exceptions to the aforementioned properties may in fact indicate suspicious or potentially deleterious mutations, which are of great interest in applied phylogenetic studies [129]. Furthermore, known or suspected mutations of this kind can be explicitly modeled as patterns and their positioning in a phylogeny assessed. In particular, those affecting important metabolic functions are expected to prevent or hinder the reproduction of the organism, and consequently should be confined in or near terminal leaves. Just as some mutations may ordinarily be forbidden altogether as global patterns, observed and feasible deleterious mutations may be permitted subject to certain restrictions. Specifically, it may be demanded that, if a hazardous pattern appears, it has no offspring, i.e., it is a leaf in the phylogeny; or, to provide some flexibility, it may be allowed that all descendants, if any, are reached in at most k steps (AXk): terminal (p)≡AG (p→leaf) (4.17) terminal (p, k)≡AG p→AXk(leaf)(4.18) In this case, leaves (self-loops in the Kripke structure) must be detected without reference to any particular sequence. This is easily achieved by performing an equality comparison between the valuations of AP of the target state and all its successors: leaf ≡Vp∈AP p↔AX (p) (4.19) This last example is representative of properties that perform conditional explorations of the phylogeny. Lineage-specific verification represents a further step forward, where patterns would be used to define the sets of states of interest. Notice that pattern-based checking of haplogroup classifications falls within this category. Single nucleotide polimorfisms (SNP’s) define genetic markers that differentiate members of the same biological specie. 4.4. A COMBINATION OF TREE AND SEQUENCE PROPERTIES 45 4.4. A Combination of Tree and Sequence Properties Some properties do not fit exclusively into one of the previous classifications but are a mix of complex sequence and tree properties. At first sight, it is reasonable to consider relatively simple properties based on a tree topology and an associated sequence alignment, exemplified by the following real properties. Suppose for now that the alignment comprises a number of cladistic characters indexed 1 through l, sequences seq are words of length lover an alphabet Σ, and seq[i] = σ(seq[i]6=σ) means that σ∈Σ appears (not) in position iin a state sequence (recall the AP definition in Definition 3). For readability and compactness in these examples, we will refer to the atomic propositions seq[i] = σ(seq[i]6=σ) of a state sequence as σi(σi). In the case of finite domains, such as the set of DNA sequences, the evaluation of logical quantifiers ∀and ∃can be substituted by multiple instances of boolean formulas connected by the Vand Woperators. Consider the following example. It determines whether a given tree is free of back mutations, which we abbreviate BM (equivalently, it detects those points in the tree where back mutations occur, if any). detectBM ≡Vl j=1AG (¬hasBM (j)) (4.20) hasBM (col)≡_ σ∈Σ σcol ∧EF (σcol ∧EF (σcol )) (4.21) In these two formulas we present a non-trivial modeling example of a cladistic property with a heavy use of sequence data. The goal is to detect back mutations in the tree, which may be encouraged by the recovery of an ancestral phenotype or a random walk in a low selection pressure zone in the DNA [65, 1]. To this end, we need to formalize the concept of back mutation: given a node in the tree (which itself defines a subtree) and an alphabet Σ, there is a back mutation in that subtree involving a position of the alignment, say j, if at some point in some descending path from the node we find in position ja different symbol (σj) than that found in the root of the subtree (σj), and if at some point in the subtree hanging from that intermediate the symbol from the root reoccurs. The formula hasBM (col) models this condition by nesting EF operators (a node satisfies a property which eventually some other descendant does not satisfy, but is fulfilled once again at some point in the future) and repeating the check for every symbol that may occur in the node. Finally, the global formula detectBM iterates the check over the positions of the alignment and extends it to all tree nodes. The second example detects a different family of complex mutations: those which affect groups of consecutive positions, a pattern of particular relevance to 46 CHAPTER 4. TEMPORAL LOGIC protein-coding regions. We dub this multiple nucleotide polymorphism (MNP), as opposed to SNP. detectMNP ≡Vl−1 j=1 AG (¬startsMNP (j)) (4.22) startsMNP (col)≡Wσ,τ,υ∈Σ(σcol−1∧τcol ∧υcol+1 ∧(4.23) EX (σcol−1∧τcol ∧υcol+1)) where for brevity the special treatment required when col = 1 has been omitted. In this case, σ,τand υare characters of the state sequence seq in positions col, col −1and col +1respectively. The property startsMNP can be extended to accept these characters as input parameters. It is clear that any verification of properties from the alignment (e.g., covariation) can be fulfilled in the phylogeny. Generally speaking, haplogroups are defined as clades whose members share a common mutation denoted by a SNP or MNP. The parametrization startsMNP (col, σ, τ, υ) will work as the membership function hidefined in Equation 4.9–4.11 as it will label the clade root (MRCA) which has the nucleotides σ,τ,υaround the position col. Given a set of clades and a known MNP, the detection of the most likely haplogroup will consist of the selection of the root with the greatest subtree satisfying the detectMNP property. Next, we must check the conservation of the MNP for all the populations inside the clade so that no spurious or external haplogroups (startsMNP) corrupt the current subtree. Using model checking packages, a fixed point algorithm will symbolically compute the set of states satisfying startsMNP for each candidate of column and characters. The objective for the haplogroup fitness is the selection of the ancestor with the maximum group of descendants verifying the property. Other related biological properties would have different maximization (minimization) purposes. Chapter 5 Tools and Experiments -“Alright,” said Deep Thought. “The Answer to the Great Question...” -“Yes...!” -“Of Life, the Universe and Everything...” said Deep Thought. -“Yes...!” -“Is...” said Deep Thought, and paused. -“Yes...!” -“Is...” -“Yes...!!!...?” -“Forty-two,” answered Deep Thought, with infinite majesty and calm. -“Forty-two!” yelled Loonquawl. “Is that all you’ve got to show for seven and a half million years’ work?” -“I checked it very thoroughly,” said the computer, “and that quite definitely is the answer. I think the problem, to be quite honest with you, is that you’ve never actually known what the question is.” —Douglas Adams, The Hitchhiker’s Guide to the Galaxy 5.1. Introduction The use of model checking techniques is domain-independent: given a phylogenetic tree and a specification of its biological properties, the verification software automatically checks the correctness of the system. In the event of failure to comply with the specification, the software outputs the scenarios which infringe the property as counterexamples. 47 48 CHAPTER 5. TOOLS AND EXPERIMENTS This chapter starts introducing a list of model checking tools adapted for different temporal logics. It also describes the implementation of the branching-time phylogenetic tree (Definition 3) into a model checking verification tool and evaluates its performance using phylogenetic properties. More in detail, we evaluate a set of back mutation properties over a real data set with the NuSMV tool [45]. In future sections and chapters we will also use the model checking software PRISM [124]. Although there exist multiple model checking tools as explained in Section 5.2, NuSMV and PRISM suffice for our requirements. 5.2. Model Checking Tools The consideration of the different types of temporal logics involves the implementation of a pool of model checker packages for supporting them (see Table 5.1). With respect to stochastic and quantitative temporal logics, PRISM is a model checker adapted for “formal modeling and analysis of systems that exhibit random or probabilistic behavior” [124]. In addition, PRISM has been successfully applied in a nearby application domain such as biological pathways [103, 123, 104, 73]. It seems suitable for the statistical analysis and quantitative evaluation of properties (e.g., mutations). UPPAAL is a commercial software package that offers “an integrated tool environment for modeling, validation and verification of real time systems modeled as networks of timed automata” [21] . This feature can be exploited in biology so as to provide an explicit definition of evolution clocks. Furthermore, NuSMV [45], a powerful tool for verifying CTL formulas, includes quantitative extensions related to the computation of paths with maximum (minimum) lengths between sets of nodes satisfying a property. It also includes a simple database manager for reusing the results of previous properties. In addition, the classic LTL open source model checker SPIN offers a framework that can be tuned as desired [93]. Sliced model checking presented in Section 7.3 can be coded over these two packages [139, 143]. In conjunction with the previous logics, on-the-fly model checking improves performance by exploring the state space on demand [24]. One of the best model checkers in this field is PROD [169]. Moreover, multicore and network-oriented tools such as DiViNe [15] (for LTL), Murphy [126] or PVeStA [5] (stochastic) are available in order to speed up efficiency. Finally, TLQSolver is a model checker that manages CTL queries and allows the mining of properties that match a specific pattern [40]. The definition of patterns helps to select sets of nodes of the phylogenetic tree satisfying an abstract relation in a similar way to an SQL query. These result sets can then be analyzed in order to extract the differences among them. 5.2. MODEL CHECKING TOOLS 49 Name Main feature Properties Language Platform NuSMV [45] OpenSource LTL, CTL, RTCTL, PSL Windows, Unix, MacOS PROD [169] On-the-fly Model Checker CTL Linux SPIN [93] Generic Model Checker LTL Windows, Unix DiViNe Tool [15] Distributed Multicore Model Checker LTL Unix Eddy Murphi [127] Distributed Multicore Model Checker Assertions Unix PVeSta [5] Statistical Parallel & Multicore Model Checker PCTL Windows, Linux, MacOS PRISM [124] Probabilistic & Quantitative Logics PCTL, CSL, LTL, PCTL* Windows, Linux, MacOS UPPAAL [21] Commercial Software for Real Time Systems TCTL Windows, Linux TLQSolver [40] Temporal Logic Query Checker for Mining Properties Query CTL Linux Table 5.1: List of some available model checkers. 50 CHAPTER 5. TOOLS AND EXPERIMENTS Figure 5.1: Boxed nodes indicate the states where the corresponding property is held. 5.3. Codification of Branching-time Phylogenies We have used two different kind of model checkers for the implementation of the branching-time phylogenetic tree. First of all, NuSMV [45], a well-known model checking software tool compatible with Cadence SMV [125, 138], is used for the evaluation of qualitative properties defined in temporal logic. Secondly, PRISM [124] is used for the evaluation of quantitative and probabilistic properties in Part IV–V. They are both freely available and they support the implementation of the whole spectrum of biological properties defined in Table 4.1. A description of the Kripke structure and the atomic propositions in the NuSMV syntax must be provided by the user as the input for the model checker. To this end, we precompute the sequence alignment and a phylogenetic tree. The translation from the phylogenetic tree to the NuSMV syntax has been performed automatically by a BioPerl script [153]. The script can be upgraded in order to include extra features such as the generation of multiple instances of phylogenetic trees, bootstrapping and so on. The inclusion of vector constructs, macros and a rich support for various logics operators in the NuSMV syntax facilitates a compact characterization of DNA and protein sequences as strings of characters. Figure 5.2 shows the implementation of the branching-time phylogenetic tree of Figure 2.1 in SMV code. The main module describes the topology of the evolutionary tree, where the names of the tree nodes (taxa) are defined symbolically 5.4. PERFORMANCE RESULTS 51 (X, . . . , S), and id variables label the states. In this case, only the DNA string is defined but extra information could be considered in future works. The init and next clauses are used to mark the root of the tree and the successors of a given state. The second part of the description consists of a function returning the DNA string associated to each node. Usually, model checking packages store the previous description of the transition system in an Ordered Reduced Binary Decision Diagram (ORBDD) [33]. These data structures are efficient graph representations for characterizing and manipulating sets or relations using boolean functions. Internally, each state of the Kripke structure is characterized by a set of atomic propositions showing the DNA associated to that state (i.e., the taxon sequence seqiσ ≡(seq[i] = σ)). Every atomic proposition can be expressed with a boolean variable indicating whether or not the character σis present in the ith position of the DNA sequence. Thus, those combinations that make true a boolean function f(seq1A, seq1C, . . . , seqlT ) will identify the DNA sequences with length lbelonging to the phylogenetic tree. This feature provides an easy scheme for the manipulation of sets in a symbolic way because a string can be identified by means of a logical expression. It should be noted that the verification of a phylogenetic formula is intimately related to the identification and manipulation of sets satisfying a boolean function characterized in temporal logic. In fact, the last step of the verification process often consists of checking that the initial state (the tree root) belongs to the result set. However, the ordering of the boolean variables in the representation of the characteristic function has a high influence in the data structure and it determines the size (number of nodes) of the diagram. In the next section, we verify phylogenetic properties using the description depicted in Figure 5.2 as input data. We aim to discover the performance trends and analyze the impact of the variable ordering. 5.4. Performance Results for Monolithic Verification of Phylogenetic Properties The performance evaluation of our system has been measured with human protein alignments retrieved from GenBank [23]. In particular, we selected genes of respiratory complex I encoded in mitochondrial DNA (mtDNA). We have chosen them because they are biologically interesting and varied in length, which makes them suitable for a complete performance analysis. This data set includes ND5, one of the biggest genes in mtDNA. Thus, the experimental results will define approximate upper bounds which can be used as a sound reference for mtDNA experiments. All tests have been run on a scientific workstation Intel Core 2 Duo E6750 @ 2.66 GHz, 8 GB RAM and Linux. Notice that the NuSMV uses a single 52 CHAPTER 5. TOOLS AND EXPERIMENTS core. We start analyzing the time and memory usage for the construction of the phylogenetic Kripke structure labeled with protein sequences. Table 5.2 and Table 5.3 show the time and memory consumption with respect to the sequence length and alignment size. Alignment Size 750 1000 1250 1500 1750 2000 ND4L (98) 13.5 27.6 46.0 67.4 89.6 124.7 ND3 (115) 7.0 31.4 53.6 78.1 105.7 142.8 ND6 (174) 18.1 26.3 71.6 123.2 160.8 205.9 ND1 (318) 8.8 13.2 84.2 203.3 246.7 371.6 ND2 (347) 10.8 13.9 75.8 217.5 322.5 420.0 ND4 (459) 14.2 22.4 29.4 191.3 417.4 499.0 ND5 (603) 19.0 92.1 106.8 314.9 518.5 728.8 Table 5.2: Seconds needed for the creation of the Kripke structure and the storage of protein sequences. Alignment Size 750 1000 1250 1500 1750 2000 ND4L (98) 102 123 146 170 194 211 ND3 (115) 111 135 161 193 218 244 ND6 (174) 161 199 233 267 300 351 ND1 (318) 275 349 409 472 533 595 ND2 (347) 307 388 454 520 586 679 ND4 (459) 423 512 598 712 798 890 ND5 (603) 568 682 823 940 1054 1168 Table 5.3: Megabytes needed for the creation of the Kripke structure and the storage of protein sequences. The time increases quadratically with the number of sequences and linearly with the gene length, except for some erratic behavior due to the use of BDD diagrams (Figure 5.3). On the other hand, the memory usage increases sublinearly with the number of sequences and the gene length, which is very encouraging from a computational point of view. It is possible that these moderate trends are partly due to the use of highly conserved genes and closely related sequences. Nevertheless, the huge amount of time and memory required for the worst case indicates that the codification of data should be optimized. Chapter 6 Conclusions Logic will get you from A to B. Imagination will take you everywhere. —Albert Einstein After the initial motivation and the introduction of the technical and conceptual basis in Part I, the next step consists of the specification and evaluation of phylogenetic properties. In the first chapter of this part, we have presented a classification of common phylogenetic properties according to the information they inspect. In particular, we have detected three kind of formulas: a) topological or tree properties, that need the structural information of the phylogeny, b) sequence or state properties, that work with the sequences attached to the nodes, or c) a combination of both. Moreover, we give a brief tutorial showing how to proceed methodologically for the incremental building and composition of specifications for each family. A taxonomy of important biological properties is summarized in Table 4.1, together with their temporal logic formulas. The inclusion of explicit time distances in certain properties motivates the addition of quantitative features in future sections. The second chapter of this part is devoted to the experimentation. We start with the state of the art of model checking tools and the selection of the best one that fits for our requirements (mainly, the support of branching-time temporal logics such as CTL). To this end, we select NuSMV although others such as PRISM are suitable for future studies in Part IV–V. We continue with the modeling of the phylogeny and the specification of properties in NuSMV syntax. During the experimentation, we have evaluated formulas that search for back mutations in the tree. We have seen that the initialization phase (creation of the associated Kripke structure) is much costlier than the verification process of a single formula. The experimental results show that the initialization time increases quadratically with 59 60 CHAPTER 6. CONCLUSIONS the alignment size and linearly with the sequence length. Additionally, memory consumption is sublinear in both cases. Despite the resources needed for model checking, it is still competitive thanks to the symbolic manipulation of huge amounts of data. The data set used there consists of proteins coded by genes from the mtDNA genome, which are substantially smaller than those from nuclear DNA. As the temporal cost increases mostly with respect to the sequence length, the phylogenetic analysis of large genes and genomes becomes the major bottleneck. Thus, scaling the model checking verification process both in time and memory is one of our research priorities for the next parts. The difficulty of implementing a completely new and adapted software, the existence of few (but good and powerful) model checking tools such as NuSMV or PRISM, and the availability of its source code favors the reuse and customization of current tools instead of concerning with new implementations. Part III Model Checking Adapted to the Phylogenetic Context: Introducing Concurrency and Solving some Bottlenecks From all the human genome, approximately the 99.9995% belongs to nuclear DNA whilst the remaining 0.0005% corresponds to the mitochondrial DNA used for the experimentation in Section 5.4. The cells of some salamanders may contain up to 40 times more DNA than those of humans [83]. Conversely, the number of available records in GenBank is currently doubling approximately every 35 months [23]. Hence, the extension and scalability of model checking for bigger alignments and phylogenies represents a dramatic question. Here, we try to solve all these problems by applying an strategy of divide and conquer. This part contains two chapters. The first one is devoted to the presentation of technical and algorithmic solutions for the problem of model checking with large genomes. This includes the use of external databases for storing the sequences and the distribution of model checking computations. Instead of saving the complete genome as atomic propositions on the phylogenetic tree, the states are labeled with pointers to the database, leading to a light memory-consuming version of the Kripke structure. Additionally, distributed model checking is divided in two subsections, explaining both the structural partitioning of the tree and the slicing of the states with DNA as feasible options for scaling the system. These solutions are compatible and they can be applied together. The structural partitioning of the tree attacks the problem by grouping the nodes in chunks and sending the partial data structure and specifications to computational units. Conversely, the slicing of the states maintains the Kripke structure but generates multiple copies of the original tree, each one with a portion of the DNA. After this, the second chapter faces the implementation of all these techniques and integrates them inside a workflow. We evaluate a set of back mutation properties over a real phylogeny with thousands of nodes, increasing the complexity with respect to previous experiments. The experimentation results show the feasibility of our approach, showing competitive time and memory requirements. 63 Chapter 7 Sliced and Distributed Model Checking Nothing is particularly hard if you divide it into small jobs. —Henry Ford 7.1. Introduction As a prominent advantage, model checking allows us to uncouple software tools from the definition of properties and it also hides the underlying implementation technology. Besides, these properties can be exported and evaluated in other structures (i.e., trees or networks) so as to compare the results and define metrics. Nevertheless, the performance is penalized when scaling the system for large phylogenies and alignments (Section 5.4). The underlying problem of standard model checking is the great amount of phylogenetic data it has to deal with: the information associated to each node of the tree is strongly related to the DNA sequence of the specie (up to millions of nucleotides). Currently, close to 2 ×106 species have been cataloged [110], but genomes exhibit an enormous range of sizes and complexity, with Homo sapiens ranking at 3 ×109base pairs. In order to cope with this problem, two strategies are presented here. The first one consists of partitioning the graph structure into a set of subgraphs, each one representing a subproblem of verification so as to speed up the computation time and distribute the memory consumption. Two subtactics are considered depending on the division method: the partition of the tree into subtrees, each one managed by a different model checker (Section 7.3); or the slicing of the tree, each slice containing a copy of the original tree but only a portion of the DNA sequence (Section 7.3). The generic techniques based on distributed model checking were 65 66 CHAPTER 7. SLICED AND DISTRIBUTED MODEL CHECKING presented in [31], but the inherent acyclicness and structure of trees facilitates the introduction of new simplifications and adaptations that has been barely developed in the community yet. We also present the novelty of sliced model checking as an adaptation to the context of phylogenetics. The second strategy is based on uncoupling the DNA information of the phylogenetic tree and exporting the alignment to an external tool specialized in the management of large information systems, for example, a database (Section 7.4). This leads to a light tree structure labeled with pointers to the elements of the database that can be efficiently manipulated by current model checking tools. The database strategy is firstly introduced here for the domain of phylogenetic model checking. In this chapter, we explain the best of the two worlds. We integrate these methodologies (distribution and databases) and show that a combination of both strategies allows us to work with real phylogenies. The introduction of all these methodologies need the presentation of the sequential model checking algorithms for the contextualization. 7.2. Model Checking Algorithms Current model checking tools load the model and the temporal logic formulas, evaluate them and compute a counterexample if so required [63, 13]. It is possible to reuse known properties and even results, if these are available. Evaluation results may discover new meaningful information that will be reused in subsequent refinements of the phylogenetic tree. Verification of temporal formulas is formalized under the framework of set theory; we follow this convention throughout the section. In fact, the traditional model checking algorithm is usually presented as a recursive function which computes the set of states satisfying a CTL formula in a Kripke structure (Algorithm 1). In order to evaluate the temporal operators EG (ψ) and E[ψ1Uψ2], the greatest and least fixed points are computed, respectively. Both fixed point sets can be obtained as the result of a breadth-first search (Algorithm 2). In particular, the call fixedpoint (M, Sat (M, ψ),∅, S) produces the greatest fixed point, and fixedpoint (M, Sat (M, ψ2), Sat (M, ψ1),∅) produces the least fixed point. The time complexity of verifying a CTL formula φagainst a Kripke structure is linear in |φ|and the size of the Kripke structure, with |S|the number of states, |R| the number of transitions and |φ|the number of logical connectives and temporal operators of the formula [63, 149]. More generally, the complexity is Θ((|S|+|R|)∗ |φ|). The decomposition of a specification for the verification of the integral parts favors the distribution of the computations in next generation algorithms. 7.3. DISTRIBUTED MODEL CHECKING 67 Algorithm 1 Algorithm Sat (M, φ) Require: M= (S, S0, R, L, AP) is a Kripke structure Require: φis a CTL formula Ensure: A subset of states of Sthat satisfies φ if φ≡ > then return S{Set of states from the Kripke structure M} else if φ≡p∈AP return {s:p∈L(s)} else if φ≡ ¬ψreturn S\Sat (M, ψ) else if φ≡ψ1∨ψ2return Sat (M, ψ1)∪Sat (M, ψ2) else if φ≡EX (ψ)return {s: (s, s0)∈R, s0∈Sat (M, ψ)} else if φ≡EG (ψ)return fixedpoint (M, Sat (M, ψ),∅, S) else if φ≡E[ψ1Uψ2]return fixedpoint (M, Sat (M, ψ2), Sat (M, ψ1),∅) end if Algorithm 2 Algorithm fixedpoint (M, Sat (M, φ), Sat (M, ψ),Init) Require: M= (S, S0, R, L, AP) is a Kripke structure Require: Sat (M, ψ) and Init are the sets of initial states (returned by calls to Sat algorithm) and Sat (M, φ) is the set of final states Ensure: A set of states that represents paths going from Init (or Sat (M, ψ)) to Sat (M, φ) New ←Init repeat Old ←New New ←Sat (M, ψ)∪(Sat (M, φ)∩ {s: (s, s0)∈R, s0∈Old}) until New =Old return New 7.3. Distributed Model Checking State of the Art Sometimes, the storage in local memory of the phylogenetic model together with the atomic propositions (both initial and new discovered properties) leads to a memory bottleneck in the model checker tool. It appears when characterizing taxa from phylogenetic trees as states in Kripke structures, that is, when labeling states of the transition system with long DNA or protein sequences (Section 5.4). In addition, the characterization of properties with temporal logics produce large sets of atomic propositions to be verified sparsely. Thus, memory usage arises as a major limiting factor in the analysis of complex systems, often in association with vast state spaces, and therefore with long execution times. In this regard, the conventional monolithic techniques conceived 68 CHAPTER 7. SLICED AND DISTRIBUTED MODEL CHECKING for improving the performance in industrial model checking applications fail to manipulate the phylogenetic tree in an efficient way due to the inherent features of the biological data, mainly the huge ratio of labels per state (i.e., DNA sequence). This problem persists in spite of the many methods which have been devised to scale model checking procedures. Of these, symbolic model checking is perhaps the most widespread [99, 85]. There exist several general-purpose memory techniques that can be applied to alleviate the problem of memory footprint, such as abstractions, partial reductions and symmetries (for a review, see e.g. [63, 31]). Current efforts in this area revolve around two main topics. Firstly, compositional reasoning [47], itself a classic approach based on the verification of local properties associated to each of a collection of “components” (e.g., codons and genes in the context of biological sequences), proceeding incrementally to infer global properties of the system through a bottom-up strategy. Chief in importance among these techniques is the assume-guarantee paradigm [135], which operates by establishing a collection of assumptions about the environment of a component and verifying the latter subject to the former. Alternative approaches are exemplified by [78], where they focus the trouble from the point of view of temporal logic formula decomposition. However, all these methods are ineffective when applied automatically in isolation to tightly coupled systems made up from highly interdependent components [48]. Only a few theoretical works have been published in relation to the complexity and advantages of compositional model checking over these simple structures [39]. Secondly, we have methods that exploit the explosive availability of multicore (shared-memory) computers [96], fast interconnection networks [165] or MapReduce in clusters [22]. Generally, these operate by partitioning the system as defined by the Kripke structure and distributing the chunks among available computing units (both storage of the partial Kripke structure and computation of satisfiability of logic formulas [84, 30]). These would be applicable here, yet by themselves they are ineffective, as they address the size of the structure in number of states and not the complexity of each state, which is the other limiting factor in phylogenetic model checking. An approximation to this last group of methods is considered in this section, but adapting them to the special context of trees. We also give the notions of sliced model checking, an innovation that focuses on the complexity of the atomic propositions for the tree states and complements the partition and distribution of the phylogenetic tree in subtrees. Alternatively, another strategy for phylogenetic model checking would consist of uncoupling atomic propositions from each state of the transition system and a) distribute them among several instances of a model checker tool (Section 7.3), or b) store the atomic propositions in an external database (Section 7.4). 7.3. DISTRIBUTED MODEL CHECKING 75 for the tree nodes having the nucleotide σin the ith position of the DNA sequence. The classic model checking algorithm parses the formula and starts verifying the inner subformulas first. In the derivation tree of the CTL grammar for the given formula, we reach a propositional operator at some point of the recursion (boxed lines in Algorithm 3). The verification of formulas with propositional operators, such as ψ1∨ψ2, begins with the computation of the satisfiability sets Sat(M, ψj), j={1,2}. These sets allow us to distribute the computation in parallel, as long as we compose the partial results with a synchronized union. The support of ψjis kψjk={pi∈APi|pior ¬piappears in ψj}. The verification of ψjis mapped to a remote model checking tool that has a copy of the phylogenetic tree labeled with the slice TI=Si∈IAPi, with I⊆ {1. . . l}a set of index, l=length(DNA) and AP =Sl i=1 APiin this case. The remote model checker also executes Algorithm 3 but with the set of atomic propositions belonging to kψjk ⊆ TIj. In order to obtain a perfect distribution, Tj={1,2}kψjk=∅. The notation par(TIj, Sat(M, ψj)) means that Sat(M, ψj) is computed in parallel in the remote model checker TIjassociated to the slice containing kψjk.LTIis the labeling function of the states of that slice. The computation of the CTL paths is denoted with the fixed point algorithms µZ/υZ. The number of parallel instances of the model checking algorithm is determined by the granularity in the slice of the AP set, but only one of those (i.e., the “coordinator”) receives the full phylogenetic property as input. Algorithm 3 Algorithm Sat (M, φ) Require: M= (S, S0, R, L, AP) is a Kripke structure Require: φ(p1, p2, . . . , pl) is a CTL formula Ensure: A subset of states of Sthat satisfies φ if φ≡ > return S{Set of states from the Kripke structure M} else if φ≡pi∈AP return par(TI,{s:pi∈LTI(s)}) with i∈I else if φ≡ ¬ψreturn S\par(TI, Sat(M, ψ)) else if φ≡ψ1∨ψ2return par(TI1, Sat(M, ψ1)) ∪par(RI2, Sat(M, ψ2)) else if φ≡EX (ψ)return {s: (s, s0)∈R, s0∈Sat (M, ψ)} else if φ≡EG (ψ)return υZ. (Sat (M, ψ)∩EX (Z)) else if φ≡E[ψ1Uψ2]return µZ.(par(TI1, Sat(M, ψ1)) ∪par(TI2, Sat(M, ψ2)) ∩EX (Z)) end if By now, we consider the access to the atomic propositions transparent to the underlying technology. However, we must advance that this is usually the most time-consuming part during the experimentations in Chapter 8. This fact moti- 76 CHAPTER 7. SLICED AND DISTRIBUTED MODEL CHECKING vates the introduction of information systems optimized for the management of huge amounts of phylogenetic data in the next section. Finally, the size of the slices depends on the target DNA regions we desire to analyze (i.e., single nucleotides, genes, chromosomes, . . .), the kind of properties we want to verify and the hardware requirements we have. A high number of slices will provide of a better level of parallelism and low hardware requirements (CPU, memory), but it will be limited by the potential appearance of bottlenecks during the composition of results. The detection of back mutations is a degenerated example of this case as it only needs a tree labeled with a single nucleotide. More intricate properties such as the detection of base correlations in a respiratory pathway advices the introduction of sliced model checking for analyzing each protein gene in isolation. 7.4. Model Checking Using Databases The second step of our approach consists of uncoupling the atomic propositions from the model checker. The use of external databases as a repository of biological sequences alleviates the memory explosion problem when the local storage of trees with partial DNA is not enough. Moreover, the database manager simplifies the interface to the DNA data because it usually allows concurrent queries and it hides the internal synchronization and data structures. In a general sense, trees are labeled with pointers to DNA sequences stored in an external server. As the AP set is now structured as consequence of Section 7.3, the atomic propositions can be purely mapped into a relational database. For example, a generic sequence seq =σ1σ2. . . σlis composed of a set of elements seq[i] = σi (in terms of AP notation) whose σivalue is univocally assigned to the associated column seq[i] in a relational table. The DNA alignment is stored in a single relational table with a row per sequence. Each row has two important fields: an identifier (i.e., the GenBank accession [23]) plus the plain string of nucleotides. In fact, the tables of relational databases can be seen as matrices that we can slice by row (number of taxa) or by column (dividing the DNA in substrings) in case of the atomic propositions were stored in separated columns. These (sub)tables would be stored or replicated in separated hard discs, servers or cluster of nodes so as to allow parallel access to the DNA data and to improve the communication bandwidth between the database and the model checker. In addition, recent versions of SQL servers support multicore CPU’s, which improves the time response when attacking the server with several queries. This approach allows to scale in memory and speeds up the system. A centralized database is source of bottlenecks during the concurrent access of multiple model checker clients: dis- 7.4. MODEL CHECKING USING DATABASES 77 tributed database managers can be used in order to optimize the performance of the transactions. Then, an external database manager is required for the storage of the AP set. In this way, the model checker algorithm will reduce the memory consumption by only saving a set of references (pointers) to the atomic propositions of each state plus a boolean variable with the truth value of the current proposition. In a hybrid system that combines slicing techniques and databases, the evaluation of a CTL formula consists of translating the last recursion call (the computation of the result set {s∈S|p∈L(s)}in Algorithm 3) into a database transaction that gets the array (s[1] = σ1, s[2] = σ2, . . . , s[l] = σl) from the appropriate record of the relational table. Each character σiof a sequence can be stored in up to a independent column of the relational table: atoms of the type (s[i] = σi) are potentially read in concurrency and the tuple is rebuilt as an “intersection” set (logical conjunction) of partial results. The projection function (Definition 6) is identifiable to the selection of a column subset. The slicing of the alignment is relegated to the vertical or horizontal partition of the database relational table. Algorithm 4 is an adaptation of Algorithm 3 where par(database, ·) is a parallel query to the database asking for the set of id’s satisfying the atomic propositions (in case of the database front-end accepts multiple queries simultaneously). It substitutes the message par(TI,·) requesting the data to the rest of slices in Algorithm 3. Thanks to the centralization of the sequences in an external repository, the parallel evaluations par(TI,·) are replaced by par(·) for the rest of situations, which can be executed in a multicore CPU using threads. Therefore, Algorithm 4 requires a shared memory machine for its execution and not a pool of sliced model checking workers. Algorithm 4 Algorithm Sat (M, φ) Require: M= (S, S0, R, L, AP) is a Kripke structure Require: φ(p1, p2, . . . , pl) is a CTL formula Ensure: A subset of states of Sthat satisfies φ if φ≡ > return S{Set of states from the Kripke structure M} else if φ≡pi∈AP return par(database, {s:pi∈LRI(s)}) with i∈I else if φ≡ ¬ψreturn S\par(Sat(M, ψ)) else if φ≡ψ1∨ψ2return par(Sat(M, ψ1)) ∪par(Sat(M, ψ2)) else if φ≡EX (ψ)return {s: (s, s0)∈R, s0∈Sat (M, ψ)} else if φ≡EG (ψ)return υZ. (Sat (M, ψ)∩EX (Z)) else if φ≡E[ψ1Uψ2]return µZ.(par(Sat(M, ψ1)) ∪par(Sat(M, ψ2)) ∩EX (Z)) end if 78 CHAPTER 7. SLICED AND DISTRIBUTED MODEL CHECKING Finally, relational databases can execute internally the model checking algorithms for the evaluation of CTL formulas. Procedural language extensions for SQL such as PL-SQL support this possibility [151]. The main advantage is that they avoid the exportation of DNA data from the database to a model checking tool and the associated bandwidth bottleneck. The lack of symbolic data structures that are habitual in model checking tools such as ORBDDs [33] is resolved by the manipulation of the data with the inherent adaptations for databases. An evaluation of SQL model checkers for phylogenetic data will comprise our future work. As far as we know, this is a novel approach because model checking was mainly used for verifying database correctness [43] and querying temporal databases [60, 44], but databases has never been used as an external repository. Furthermore, relational databases cannot only store phylogenetic data, but also partial verification results or counterexamples. Chapter 8 Workflow The White Rabbit put on his spectacles. “Where shall I begin, please your Majesty?” he asked. “Begin at the beginning,” the King said gravely, “and go on till you come to the end: then stop.” —Lewis Carroll, Alice in Wonderland 8.1. Introduction In this chapter, we describe the implementation of the different methods presented previously. We compose them in a workflow by integrating the modules encapsulating each technique. Here, we detail the internal characteristics of our framework as well as the performance results with large phylogenetic data. The combination of all these optimizations facilitate the efficient management of big phylogenies by current model checking tools. In particular, we develop an hybrid system: we code a multi-threading version of Algorithm 4 and use an external database for saving the sequences. To this end, we modify the model checker NuSMV and send queries to a MySQL server. This part is thus divided in three points. Firstly, Section 8.2 draws a panoramic overview of the pipeline, details every component and how they are connected through its interface. Next, Section 8.3 gives brief notions about the experimental results and performance of model checking with these new approximations. We search for the presence of back mutations over ZARAMIT [27] (a real phylogenetic tree) and evaluate the temporal costs. Finally, Chapter 9 summarizes the conclusions of this part and suggests the future work. 79 80 CHAPTER 8. WORKFLOW 8.2. Description Figure 8.1 represents a graphical description of our workflow. The central core of our approach is the model checker package NuSMV v2.5.4 [45], which is a well-known public software. It is surrounded by a set of modules and tools that accommodate the phylogenetic data for each step of the verification. The system also counts with an external MySQL server v5.5.24-7 for saving the atomic propositions (i.e., the DNA alignment). The use of databases simplifies the coding of the sliced model checking algorithms. Here, we present an hybrid system where the new sliced model checking algorithms are highly coupled with the database (Algorithm 4). The sliced model checking algorithms are implemented in two steps. First of all, an external script analyzes the specification of the phylogenetic properties that the system must evaluate, and subsequently accesses in parallel to the atomic propositions (DNA sequences) stored in the database. After the extraction of the GenBank identifiers satisfying those AP, these sets are converted into the internal data representation of NuSMV (i.e., ORBDDs). Finally, a multi-threading version of Algorithm 4 collects the results and compounds the CTL paths for evaluating the truth of the property. In sum, our framework is thus divided in three important modules. First of all, the loader consists of a BioPython script that creates and initializes the relational database with the DNA sequences. It is executed only once during the initialization phase for loading the sequences into the database. The database habitually remains constant for the rest of the verification process, but it can be updated periodically with the addition of new aligned strings. By default, we use BioSQL, a shared database schema for storing sequences in SQL servers that is supported by several script languages [121]. Secondly, the property transformer is a BioPython script that pre-processes the phylogenetic formulas and rewrite them in terms of GenBank identifiers. Every time the user asks for the states satisfying a pattern seq[i] = σin the DNA, the script retrieves the set of identifiers of those species that satisfy the formula in the database. This script translates the phylogenetic properties expressed in terms of sequences and nucleotides into phylogenetic properties written in terms of taxon identifiers. For example, the following formula reinterprets the back mutation property (Equation 4.20): BM ≡(id ∈seqiσ)∧EF [id /∈seqiσ ∧EF (id ∈seqiσ)] (8.1) where id is the identifier of the state reached by each temporal operator EF (at the beginning, id corresponds to the identifier of the initial state). The set seqiσ ={id1, . . . , idn}contains the identifiers of the set of states satisfying the 8.2. DESCRIPTION 81 property seq[i] = σin the external database. The states of the phylogenetic tree are labeled with these identifiers. The script is optimized for sending concurrent queries to the database in order to take advantage of current multicore servers, each question demanding the set of identifiers satisfying seq[i] = σ. As result, it outputs an specification file with the original phylogenetic properties expressed in terms of GenBank identifiers. In addition, it generates a set of auxiliary files corresponding to the projections of the formula in the sliced model checking methodology. These auxiliary files, one per seq[i] = σ, store the set of identifiers seqiσ satisfying an atomic proposition of the original formula. Next, the NuSMV transformer is a BioPerl script that translates the phylogenetic tree into the NuSMV syntax (Section 5.3) and labels the states with the associated identifiers of the database. The script is easily extended in order to divide automatically the original tree in multiple subtrees as a part of the implementation of partitioned model checking, but the transformation of the equations together with the composition of counterexamples and results has to be realized manually at the moment. The subtrees can be stored either in separated NuSMV model files in order to evaluate the specifications in concurrency; or in individual modules inside the same file for a sequential evaluation and composition of the results. The cost in time/memory increases fast with the size of the tree and therefore a sequential evaluation of the specifications in the subtrees also becomes profitable in terms of improving the efficiency. The CTL parser in NuSMV may assist to the unfolding of the temporal logic formulas according to the explanation in Section 7.3. Once the phylogenetic tree is translated into the model checker syntax and the properties are translated and projected according to the slicing criterion, the workflow starts the verification process. To this end, a set of independent instances of NuSMV are launched. Each one works with a copy of the Kripke structure and a DNA projection, that is, the auxiliary files generated by the property transformer. Their mission consists of converting the set of identifiers of seq[i] = σstored in the auxiliary files into the internal notation of NuSMV (an ORBDD). At the end of the process, a multi-threading version of the model checker collects the ORBDDs for the atomic propositions and reconstructs the CTL paths. The final speed up depends on the distribution and the number of formulas and slices. The results of the verification process offer an important feedback for the refinement of the original tree and the biological assumptions. 82 CHAPTER 8. WORKFLOW Figure 8.1: Workflow diagram with the alignment, phylogenetic tree and properties as input. 8.3. Experimental Results The data set of our experimentation comprises the ZARAMIT tree [27], which has been reconstructed from 7390 aligned sequences of Human mitochondrial DNA with 16,569 base pairs. In total, the number of internal plus terminal nodes is 14512. Thanks to the integration with databases, the memory consumption of the model checker is reduced: the representation of the ZARAMIT tree consumes around 50 MB. Besides a light consumption of memory requirements, the use of references to an external database reduces the initialization time of NuSMV. Figure 7.2 shows that our architecture outperforms previous works in [138]: we spend less time and memory resources for the initialization of the model checking tool with the ZARAMIT tree than for the initialization of a protein tree with 2000 tips and sequences of 500 aminoacids. Taking advantage of these optimizations, we have employed the ZARAMIT tree for the evaluation of 25 back mutations formulas using our framework. From the point of view of performance, the database is the most important point of our workflow because it connects the model checker with the DNA sequences. The database server runs in a desktop workstation (AMD Opteron @3GHz, 4GB RAM, Debian Linux). Due to the alignment size (around 350 MB for 7390 leafs plus ancestors) and the hardware available, we work with a single 8.3. EXPERIMENTAL RESULTS 83 instance of MySQL server. The database loader spends less than 2 minutes for inserting the alignment in the database at the initialization point, but it remains constant for all the verifications. The retrieval of atomic propositions from the database lasts 7.2s for 25 concurrent SQL queries (44.6s if we serialize them). Each SQL query asks for the set of identifiers associated to the DNA in each one of the 25 back mutation properties. A second desktop workstation (Intel Core 2 Duo E6750 @ 2.66GHz, 8GB RAM, Debian Linux) is devoted to the execution of the rest of the workflow. This comprises the translation of the phylogeny into the NuSMV syntax (NuSMV transformer), the conversion of the GenBank identifiers to the internal data structure of NuSMV (auxiliary files of the property transformer) and the execution of the multi-threading algorithm that accomplishes the evaluation. By now, the back mutation properties are evaluated over the whole phylogeny. The application of partitioned Kripke structures (i.e., the verification of the back mutation property independently for each direct subtree of the phylogeny) will divide the time by more than a factor of 2 for binary trees. However, we have detected that the internal representation of the Kripke structure in NuSMV penalizes the access and manipulation of big sets during the model checking process. Given a set of GenBank identifiers whose DNA sequences satisfy the atomic propositions seq[i] = σ, the time required for the extraction of the tree states with those id and their exportation to a temporal ORBDD file raises up to 1m20s for each auxiliary file of the property transformer. Nevertheless, this latency can be alleviated because the exportation of every set of GenBank identifiers to a ORBDD file is executed in parallel. Conversely, the integration of partial results (read of a ORBDD file) and the computation of CTL paths for the verification of a back mutation property only lasts 1m 2s for the multi-threading algorithm. All these temporal costs embeds the initialization time of NuSMV for the ZARAMIT tree (around 57s). In sum, the initialization costs of NuSMV raises up to 55s-1m for phylogenetic trees with thousands of species and nodes labeled with references to the database. Besides, the framework offers a competent performance for evaluating phylogenetic specifications, such as the detection of back mutations. The first step, the extraction of big sets of states satisfying a simple atomic proposition, lasts around 1m27s (SQL query + exportation to a ORBDD). The extraction of multiple atomic proposition can be effectuated in parallel, with the number of cores the only limitation to this part. Later, the integration of the ORBDD results and the computation of the CTL paths takes 1m2s, leading to 2m30s for a complete verification of a single back mutation property. The parallel evaluation of the back mutation property in the direct subtrees of ZARAMIT will reduce the verification time by more than a half. Longer unexplored phylogenetic properties can take advantage of multi-threading 84 CHAPTER 8. WORKFLOW as the CTL verification cost depends of the formula length [13]. 10.1. INTRODUCTION 91 Figure 10.1: Phylogenetic tree for the Hominoidea. The incorporation of CTL path operators with temporal bounds allows the localization of ancestors whose descendants are placed at a certain distance and are labeled with one of the mentioned tags (see Section 11.2 for a more detailed example and notation). Traditionally, the phylogenetists solve the problem of computing distances with an ad hoc script that adds the weights of a path in a enriched tree with time. The insertion of explicit time in the CTL path operators avoids the codification of these specific scripts and maintains all the advantages of model checking exposed in this thesis. The new model checking algorithms will embed the extraction and addition of distances through the traversing of the states in the phylogeny. Besides, the algorithms of qualitative model checking are unable to manipulate quantitative information in the edges of a transition system or generate an integer as output result. This is of remarkable interest when we try to extract temporal information from the tree (e.g., distance(pongo, homo)) instead of testing with particular values. In general, the new model checking algorithms should focus on a) the evaluation of properties whose truth value depends on these quantitative relations; or b) the evaluation of properties whose result is a numerical value computed using these annotations. For the first case, the descriptive aptitudes of the logics are extended to handle with the numerical information in the equations. The existing temporal logics are recycled for this field in order to cover the maximum number of phylogenetic properties querying about quantitative data. For the second case, the model checking process is modified with respect to the typical verification of a temporal logic formula: now, it will return a number instead of a boolean value. The model checking algorithms are extended with the capability for evaluating arithmetic functions over the set of states characterized by a formula of temporal logic. As explained in Section 2.4, the set of tree states are selected according to the relations imposed by a temporal logic formula, and they can be considered as objects or entities in 92 CHAPTER 10. TIME, PROBABILITIES AND QUANTIFICATION a phylogenetic tree. The studies can be enlarged with the introduction of parameters in the equations: there exist automatic tools that infer the combination of propositions and values (quantitatives or not) that make the equations true. The aim of this chapter is to review the group of non-qualitative phylogenetic properties and classify them according to a) the kind of information that is tagged in the tree/specifications for the verification process, and b) the kind of answer the verification will return (boolean or numerical) as well as the procedure it uses for its computation. The introduction of quantitative information in phylogenetic specifications will be motivated by real examples. Looking at these examples, we will identify the requisites that the future logics should have. The design of such logics will also pay attention in the main aspects of decidability: many properties expressed in linear or branching-time logics based on Presburger formulas are undecidable. Finally, the last important factor to take into account is the availability of predefined logics and the associated verification tools covering the proposals presented here. We will present the current developments of quantitative extensions in the field of model checking, apply them to the definition of new properties over phylogenies and provide to the phylogenetists a complementary framework where they could inspect numerical relations. The enlargement of the semantics of Kripke structures and the expressiveness of the temporal logics contribute to the definition and verification of such properties. The expression and evaluation of quantitative properties involving phylogenetic hypothesis is the main advantage of this chapter. 10.2. Towards a Classification of Quantitative Properties Arising in Phylogeny Evolutionary pressure of natural selection affects differently to each taxon because they are fixed to incomparable environments and ecosystems. Even, genes concerning disparate traits of the same species mutate with apparently independent speed, generating different topologies of gene trees. Moreover, current phylogenetic inference methods working with the same input not necessarily obtain equal output phylogenies [71, 176]. Consequently, evolution is not an homogeneous process, which is reflected in the variability of regions in genomic alignments and, indirectly, the trees built from these data [131, 132, 119]. Phylogenies express these biases through the tree diversity and unbalanced trees (binary or not). Therefore, the branches that are supported by the majority of the inferred phylogenies determine the main tree skeleton and the core of the evolution. The inclusion of probabilistic information in the branches of the tree appears naturally during the statistical analysis of the phylogeny: probabilities help to the detection 10.2. CLASSIFICATION OF QUANTITATIVE PROPERTIES 93 of the consensus tree [71], i.e., expressing the rate of trees confirming the topology. The computation of these probabilities in the branches can be carried out by the model checking tools using the equations of maximum likelihood estimations (MLE) [69]. To this end, we need to extend the capabilities of the software tools for evaluating quantitative functions over the set of states selected by the traditional model checking algorithms: now, they will return a number as output result instead of a boolean value. Otherwise, likelihoods also work for inspecting the probability of reaching a given state of the phylogeny, for example, the probability of reaching a population with a chronic disease from the root (i.e., the lactose intolerance [160, 92]). For this case, we introduce explicitly the probability inside the specification of the property that we are going to evaluate: the result will be a boolean telling whether the number of states satisfying the temporal logic formula is over a certain ratio or not. Commonly, it is assumed that every descendant is reached with equal probability from the root (i.e., every transition is equiprobable). Additionally, the edges are occasionally enriched with time labels or lengths for marking the differences in the speciation speed that appear in the nature. A tree with integer weights in the branches is generated as output of a distancebased phylogenetic reconstruction method [148], which builds the tree measuring the number of mismatches between each pair of genomes. Therefore, the inclusion of temporal labels in a phylogenetic tree inherits the concept of character distances from the tree inference methods. In fact, the integer tags representing the time in the edges have a dual meaning: they can be interpreted as the rate of nucleotide variation between two consecutive genomes of the tree, or the chronological time of the evolution. Besides, the knowledge of time is useful for learning complex properties about the evolution and speciation speed, for instance, the estimation of the temporal point of divergence between species [17], the diaspora of human populations [76, 36, 101, 166] or army invasions [179]. Hence, a phylogenetic tree can be labeled with many types of tags, quantitative or not: phenotypes such as chronic diseases or endemic features, geolocation of populations and migrations, and so on. These data are specific for each phylogenetic study, but we search for a generic criteria for classifying the addition of new information. Time and probability are the most common quantitative information and they are shared by almost every phylogenetic study because of the aforementioned arguments. Therefore, the modifications of the logic, data structures and algorithms will be driven by these types of data. In sum, we can organize the tags as integers (chronological or character distances) or non-integers (probabilities) for this first approach. Apart from the new annotations of the phylogeny, the quantitative information can also be placed in the result of the model checking evaluation: instead of returning a boolean, the output is an integer or real number. 94 CHAPTER 10. TIME, PROBABILITIES AND QUANTIFICATION For compatibility with the annotations of the tree, the quantitative output of the model checking algorithms can follow the classification of integers (chronological or character distances) or non-integers (probabilities) as well. 10.3. Extending the Labels of the Phylogenetic Tree and the Results of Model Checking So, the inclusion of quantitative information in the phylogenetic analysis can be boarded at two different levels: either in the definition of properties and the annotations of the tree, or in the result returned by the model checking process. If we focus on the first level, the numerical data inserted into the tree and specifications can be divided in a) timed or distance information, b) probabilistic information and c) raw quantitative information. The first two ones are mainly related to the labeling of the tree branches, while the last one is more general and involves numerical values and comparisons in the atomic propositions [137]. Take the following example. A phylogenetic tree with advanced annotations is illustrated by Figure 11.1. It represents a phylogeny whose branches are scaled proportionally to the time. The subtrees are labeled with probabilities representing either the statistical support for that clade, or the likelihood of taking that edge during the speciation. For instance, imagine that the tips represent individuals from different ethnic groups and the leaves are labeled with their DNA, an identifier and a boolean indicating if they suffer from a specific illness. Suppose that X and Y suffer from lactose intolerance, an endemic disease that causes the inability to digest milk. The rest of members are healthy. Timed information. The branch lengths are proportional to the mutation clock pointed by the segment in the legend. Each transition implicates a change of 0.05 nucleotides (5%) per site and per unit length with respect to the direct ancestor, and thus, it provides a visual mode for studying the total rate of variation of the genomic sequences along a segment. In this example, we can inspect the number of changes between the ancestor of X and Y, and the rest of healthy individuals: it indirectly tells the estimated date of appearance of the (in)tolerance. A phylogenetic tree whose branches are proportional to the number of nucleotide changes is called phylogram, and chronogram if their edges are directly proportional to the expected evolutionary time. For distinction, the basic phylogenetic trees whose branches are not scaled with respect to neither time nor character changes are simply called cladograms. Probabilistic information. The red values attached to the internal nodes of Figure 11.1 indicate the statistical consistency or support of a clade. In brief, those numbers express the percentage of inferred or bootstrapped phylogenies 10.3. EXTENDING THE LABELS OF THE PHYLOGENETIC TREE 95 whose leaves are arranged forming the same subtree. The ordering of the tips is irrelevant for the identification of subtrees [18]. An internal node labeled with 1.0 means that the descending taxa are clustered together in every inferred phylogeny constructed using the relationships imposed by the sequences. It is possible to compare phylogenies through the computation of maximum likelihoods estimations [69]. Not only do the probabilities play an important role for the comparison of trees through the evaluation of the tree coherence, but also they are useful for showing the prevalence of chronic diseases inside genetic populations. The blue values indicate the probability of moving from one node to another following that branch: by default, we assume that every node is equiprobable and the summation of probabilities for each descendant must be 1.0. Using the blue values in this example, the probability of reaching an individual with lactose intolerance from the root is 0.25 (the likelihood of the X-Y subtree is 0.5×0.5).  Raw quantitative information. The incorporation of raw quantitative information to the atomic propositions of the tree states increases the expressive power of the system. This quantitative information is essential for describing questions about biochemical properties related to the sequences, for example, asking the polarity, hydrophobia and size of the amino acids in the synthesized protein. The necessity of including these data was already pointed in Section 4.3 during the investigation of the conservation and covariation of sequences and their attributes. The addition, manipulation and comparison of numerical labels is already supported in the verification process of some current model checking tools, but it is important to take in mind the complexity of introducing Presburger and Peano arithmetics. We classify the states with raw quantitative information in a separate category in order to emphasize the distinction of the labels (genomic characters vs numerical atomic propositions).  Apart from the previous concerns that classify the numerical labels of the tree in three families (timed, probabilistic or raw quantitative information), the need for the extraction of quantitative results in some phylogenetic analysis is clear: recall the calculation of distances between populations (pongos and homo) instead of testing particular values. If we focus on the answer returned by the verification process rather than the type of tags, the second big classification of phylogenetic properties is separated in: Boolean properties, where the result is true when the property is satisfied, or counterexamples when the property is false; Numerical properties, where the result is the computation of a number according to a function applied over a group of elements, and whose most prominent exponents are the calculation of (continuous or discrete) time distances or probabilities; and, finally, 96 CHAPTER 10. TIME, PROBABILITIES AND QUANTIFICATION Figure 10.2: Phylogenetic tree labeled with quantitative information. Parametric properties, where the result is an algebraic equation, occasionally enriched with the set of integers or characters whose values make true the temporal logic formula. This equation may involve quantitative relations or genomic sequences. The study of boolean properties has been already detailed in Chapter 4. Given an initial state, they compute a set of reachable nodes fulfilling a certain feature. The application of quantitative functions to inspect the numerical relations over that set is also important. This quantification is crucial for the computation of distances (timed or not) between sets of taxa or clades. The distances between elements are mainly employed for the calculation of topological measures (e.g., tree depth and balance [130]) and phylogenetic comparative methods (e.g., RobinsonFoulds [146]). The inclusion of distances in the topological properties presented in Section 4.2 helps to study, compare and refine the tree structure [158]. These supplementary metrics are of particular interest in the case of several phylogenies sharing identical scores under the same inference method [159]. Conversely, the evaluation of numerical properties may return probabilities too; for instance, during the calculation of the statistical support of a tree [69, 71]. We can analyze the fitness of a set of trees according to all these informations. Finally, parametric model checking allows the placement of parameters either in the specifications or in the transitions of the model. The use of parameters in the specifications helps to the identification of potential valuations that satisfy a property instead of the manual inspection of the whole state space [38]. The evaluation of a parametric formula will return a set of instances that confirm the equation. Some properties such as the detection of back mutations, or the appear- 10.4. QUANTITATIVE PROPERTIES IN PHYLOGENETICS 97 ance of conserved and correlated regions of the DNA, can be expressed in terms of parameters instead of particular instances of nucleotide characters. Properties involving the verification of numerical restrictions over time, probabilities or raw quantitative data are also susceptible to the definition of upper and lower parametric bounds in the CTL path operators [32]. In case of the placement of numeric parameters in the model, the result of the evaluation is a mathematical equation representing the relations that must be fulfilled by the values. 10.4. Quantitative Properties in Phylogenetics One of the main results of this chapter is the sorting of the phylogenetic properties with quantitative requirements. Following the catalog proposed in this section, the phylogenetic properties are organized in a Cartesian product presented in Table 10.1. In the first column, we are centered in properties related with the genome and its arrangement in the tree. If we focus on properties with a boolean result, we obtain Table 4.1. In Chapter 4 we detailed the methodology for defining and evaluating properties of this genre. The use of parameters allows the extraction of the columns (i,j) and nucleotides (σi,σj) making true the equations. By definition, the output result of a qualitative formula cannot be a number. Secondly, the inclusion of explicit time admits the specification of restrictions such as distance(ape, homo)>10, where we test if the accumulated weight for the transitions between these two taxa is over or below 10. As explained previously, the verification or computation of distances between sets of populations can be carried out by extended CTL-like formulas that are evaluated over a enriched phylogeny with temporal information. These new CTL path operators collect the weight of each branch and test if the accumulated length is over the indicated value or not (Chapter 11). Conversely, the algorithms of quantitative model checking are updated for calculating integers (distances) and float numbers as output. The calculation of distances between symbolic objects is detailed in Chapter 12. Either the distance or the symbolic identifier of species in the specifications can be replaced by a parameter. Thirdly, we test if the number of states validating a property in the phylogeny is over or below a certain probabilistic bound. Several types of phylogenetic properties accept the inclusion of probabilities in their definition. For instance, we can evaluate sequence properties (e.g., checking the ratio of DNA mutations that activate the lactose tolerance in the phylogeny) or tree properties (e.g., checking the likelihood of generating a particular topology with a DNA substitution model). Similarly to the calculation of distances, the CTL path operators that are adapted for probabilities collect the likelihood of each branch and test if the accumulated probability is over the indicated value or not. Although probabilities are essentials 98 CHAPTER 10. TIME, PROBABILITIES AND QUANTIFICATION in both types of properties, the calculation of likelihoods differs in each situation because they use different models: the evaluation of probabilities over the alignments requires a phylogeny (model) labeled with likelihoods (Chapter 15), while the evaluation of the tree support for a certain topology requires a DNA substitution model (Chapter 14). Some model checking tools such as PRISM [124] have algorithms that return these float numbers as output of the evaluation. The parameters introduced in the genome for a sequence property are identical to the parameters introduced in the first column. The DNA mutation model can be also left undetermined: in this way, the model checking tool will search for the best parameter of the DNA model that generates the topology of the tree, but this field is still unexplored and needs a more detailed study. The model checking algorithms can be updated for returning the maximum or minimum probability of satisfying a property. Finally, the properties with raw quantitative information include arithmetic comparisons in the atomic propositions. The last column is equal to the first column in terms of restrictions and expressiveness, with the addition of numbers as the only novelty in the inequalities and parameters. The type of properties explained in Section 4.3 (conservation or covariation of biochemical patterns) falls within this category. The incorporation of quantitative information may increase the complexity of the associated logics and model checking algorithms. 10.5. Conclusions The main contribution of this chapter is the identification of the limitations of phylogenetic model checking for inspecting more complex properties. To this end, we have detected new quantitative requirements for the analysis of phylogenies. The phylogenetic properties are classified and organized according to the quantitative information they use: principally the time (distance) between states and probabilities. There exists many more quantitative tags, but these two are the most representative. The evaluation of hypothesis with this kind of data requires the extension of the phylogenetic tree with temporal and probabilistic annotations, a logic with a renovated syntax, and an algorithm considering this knowledge. The inclusion of quantitative information is not exclusively restricted to the definition and verification of specifications with explicit time or probabilities: now, the model checking process may return numerical results as output as well. This last approximation demands the update of the traditional verification algorithms for managing and computing specific functions over a sets of states characterized by a temporal logic formula. In sum, the next chapters explain the incorporation of time and probabilities to model checking with intuitive examples. Our goal is the application of existing 10.5. CONCLUSIONS 99 Qualitative Time Probability Raw B Table 4.1 Phylogenetic properties enriched with distances: distance(pongo, homo)>10? Probabilities over alignments: Pr(conservation(2,3,AA)) >0.8? Likelihood of a tree topology given a DNA mutation model: Pr(Jukes Cantor, topology) ≤0.4? Biochemical properties: EG((size > 110) ∧(size < 130))? N — (Im)balance metrics Max/min distances between clades, roots and leaves Length of right and left paths Symbolic distances (Robinson-Foulds) Max/min probabilities over alignments Max/min likelihood of a tree topology — P i,j∈N,σi, σj∈Σ, seq ∈Σ(j−i): covariation(i,j, σi, σj) conservation(i,j,seq) hasBM(i, σi) Phylogenetic properties enriched with distances: X∈N: distance(pongo, homo)> X A, B ∈species : distance(A, B)>10 Probabilities over alignments: i,j∈N,seq ∈Σ(j−i): Pr(conservation(i, j, seq)) > 0.8 Inferring the parameters of a DNA mutation model for a tree topology: Pr(model, topology)≤0.4 Biochemical properties X, Y ∈N: EG((size > X)∧(size < Y )) Table 10.1: Summary of the most important phylogenetic properties (Type B: boolean; N: numeric; P: parametric). 100 CHAPTER 10. TIME, PROBABILITIES AND QUANTIFICATION developments in quantitative model checking to the verification of phylogenetic questions. The structure is arranged as indicated. First of all, Chapter 11 introduces real (continuous) time in the specifications and transition systems. It presents the syntax and semantics of Timed CTL (TCTL), with an application to compute temporal distances between objects characterized symbolically in Chapter 12. We continue with the introduction of probabilities that depend either on discrete (Chapter 14) or continuous time (Chapter 15). To this end, we use previous investigations in Markov chains and logics such as Probatilisic CTL (PCTL) of Continuous Stochastic Logic (CSLTA). Finally the application of parametric model checking is presented as a future work and we finish with the conclusions of this thesis. 11.2. TIMED LOGIC AND STRUCTURE 107 •Adelay transition (hsi, νii,hsi, νi+1i). Only the clocks are updated with νi+1(x) = νi(x) + ti+1 −ti, with ti+1 > ti, and ti, ti+1 ∈R. The clock valuation νi+1 may not satisfy any firing condition τ(ri)yet. L0(hs, νi) = L(s)∪ {g∈CC(C)|νg}is the labeling function, and AP0=AP ∪CC(C)the set of atomic propositions. Now, the notation of the infinite paths πis enlarged for considering the potential infinite set of states S0(which is caused by the Cartesian product with the uncountable number of valuations for the real time clocks Eval(C)). Hence, following the discrete transitions of the transition system, πevolves with the cadence π=hs0, ν0ihs1, ν1ihs2, ν2i. . .. In this case, the state hsi, νiichanges to hsj, νi+1iwhen νisatisfies the enabling condition τ, which indirectly determines adiscretization of the dense-time paths in the timed automaton. Otherwise, the path πremains unaltered in the same state siwhile it executes a delay transition, i.e., π=hs0, ν0ihs0, ν0+t1ihs0, ν0+t2i. . . with t2> t1. The function ρ(π, t) is updated for considering the new type of π, i.e., ρ(π, 0) = hs0, ν0i. The translation of a branching-time phylogeny to that transition system is addressed by the extension of Definition 3 with the inclusion of clock restrictions τin the set of edges Eof the tree. In a phylogeny, these constraints define upper or lower temporal bounds in the connections between species, probably marking the interval of time when those taxa diverged. The self-loops in terminal nodes of the transition relation R0lead to perpetual siphons (by default, these self-loops consume 0 ticks of time). Every tree node can be also enlarged with local clock variables for monitoring the time spent until the next state. Normally, a single clock representing the chronological time through all the tree will be enough for our purposes in phylogeny, but we keep the notation of a set of clocks Cfor generality. Finally, let’s remember the examples presented in Section 11.1 that motivated the introduction of temporal logics with explicit time. For instance, we want to discover the date of divergence of the Tibetan ethnic group with respect to its closest relative population. In particular, we aim to detect when the genetic adaptation that favors the habitability in high altitude zones appeared (it marks the tribe genome as well). To this end, we define a phylogenetic property that searches for the most recent common ancestor of the Tibetans. The following TCTL Equation 11.2 with time extensions introduces the notation that we have presented here. This formula investigates, in the phylogenetic tree of Figure 11.1, if there exists an internal node (EF) whose terminal leaves are reachable in a distance less than 3 kiloyears (AF≤3ky) and the individuals they represent possess the mutations affecting to this new phenotype: EF AF≤3ky(seq[i] = j) (11.2) 108 CHAPTER 11. TIMED TRANSITION SYSTEMS AND LOGICS Figure 11.1: Phylogenetic tree and its transition system labeled with time intervals in the branches and taxon identifiers in the nodes (the DNA sequences are omitted for readability). This equation can be exported to other examples such as the study of time in the woolly mammoth lineage. Hence, TCTL is suitable for the definition of phylogenetic properties enriched with explicit time distance. The existence of software tools and model checking algorithms for managing TCTL queries provides a complete framework for the verification of biological hypothesis. The model checking algorithms needed for handling TCTL formulas are presented in the next section. 11.3. Algorithm for Timed Model Checking The algorithms of model checking are enlarged for taking into consideration the new temporal information of the transition system. Due to the infinite computational tree, the model checking process has to do a few simplifications in order to cope with this complexity. These simplifications are founded on the fact that certain states of the transition system are indistinguishable by this logic if we focus exclusively on the time, for example, when the time elapsed in two different states is the same. The underlying idea is the identification of equivalence classes according to compatible temporal constraints. To this end, we define an equivalence of clocks [7]. Definition 11 (Equivalence of Clock Assignments).Two clock valuations ν, ν0∈ Eval(C)are equivalent (ν≡ν0) if they satisfy the following conditions. 11.3. ALGORITHM FOR TIMED MODEL CHECKING 109 For each x∈C,bν(x)c=bν0(x)cor both νand ν0are greater than cx, with cxthe largest constant of xfound in the firing conditions τand b·c the integer part of the real number. For each pair of clocks x, y ∈C, such that ν(x)≤cxand ν(y)≤cy, the fractional parts (fract) of the clocks should: •fract(ν(x)) ≤fract(ν(y)) iff fract(ν0(x)) ≤fract(ν0(y)), and •fract(ν(x)) = 0 iff fract(ν0(x)) = 0. The application of this equivalency defines a finite number of clock regions in the original transition state space. In other words, it groups the infinite set of states S0into a finite set of clusters, each one fulfilling compatible time restrictions: the clock valuations ν∈Eval(C) are organized in equivalent regions (denoted as Eval(C)/≡). The equivalence class of νis denoted by [ν], i.e., the set of ν0that are equivalent to ν. As a particular case, an end class is the equivalence class satisfying x>cxfor all clocks x. An equivalence class αis a boundary class if for any ν∈αand t∈R>0,νand ν+tare not equivalent. We must define a successor function over the clock equivalences in order to capture the transitions between clock regions [7]. Definition 12 (Successor Region).Given a set of clock regions defined over Eval(C)/≡according to the equivalence of clock assignments (Definition 11), the successor of a clock region αis β(succ(α) = β) iff for each clock ν∈αthere exists t∈R>0such that ν+t∈βand ν+t0∈α∪βfor all t0< t. The state space is reordered in terms of regions defined over the transition system, classifying the states hs, νiin a two dimensional space by proposition φ and valuation ν. The state region of s0=hs, νi, denoted [s0] with s0∈S0, is defined by [s0] = hs, [ν]i={hs, ν0i | ν0∈[ν]}. Any movement among nodes inside this set doesn’t change the validity of the state proposition, which is only affected by an inter-region transition. Hence, the clock valuations Eval(C) and the state labels determine a finite number of equivalence classes that are the basis of the model checking process. The aim of the algorithm presented in this section is the execution of the classic CTL model checking algorithm over the regions of this last finite structure, called Region Transition System (RTS) and that we will define later. The objective of the following transformations is to translate the problem of solving a TCTL equation to the standard problem of solving a CTL formula. The process of verifying a formula φis divided into five big steps, which can be summarized in pseudo code as explained in Algorithm 5 [13]. In brief, it starts transforming a TCTL formula φinto a CTL equation φ0by deleting the timed 110 CHAPTER 11. TIMED TRANSITION SYSTEMS AND LOGICS Algorithm 5 Algorithm Sat (M, φ) Require: TA = (S, S0, R, L, AP, C, λ, τ) is a timed automaton Require: φis a TCTL formula Ensure: TA φ Transformation of φinto φ0by the elimination of timed restrictions ∼cof the formula, for example, in ψ1U∼cψ2 Migration of the timed restrictions ∼cto a state condition with zan external clock, for example, in φ0=ψ1U(ψ2∧(z∼c)) Determination of the equivalence classes under ≡∗ Construction of the region transition system M=RTS(TA,φ) equivalent to the timed automaton TA Application of the CTL model checking algorithm Sat (M, φ0) to check M, hs, [ν]iφ0 restrictions ∼c. The temporal information must be reincorporated to the verification process in the second step. To this end, we introduce an extra clock z /∈Cto the model and specifications for tracking the time elapsed in traversing the regions of the region transition system. The clock zis initialized to 0, it is never reset along the path and it is updated consistently with the other clocks. Next, the formulas are rewritten in terms of state conditions with zinstead of temporal restrictions over the path operators. The states of the RTS are labeled with atomic propositions representing the clock constraints that are fulfilled in that state. Later, we calculate the equivalence classes according to ≡∗. The equivalence class ≡∗is an extension of ≡that includes the set of clocks C∗=C∪ {z}. Finally, we construct the region transition system associated to the timed automaton TA and formula φ, and we apply the algorithms of CTL. A region transition System RTS(TA,φ) is a transition system whose nodes are clustered by time and by property according to the restrictions they satisfy. It depends implicitly of the formula φthat we try to verify: its maximal constants ∼care of relevance to the clock equivalence only. In case the region transition system does not depend on (the maximal constants occurring in) φ, we simply write RTS(TA). The RTS is defined as: Definition 13 (Region Transition System).Given a TCTL formula φ, the region transition system RTS = (S∗, S∗ 0, R∗, L∗, AP∗)is an enhanced transition system TS = (S0, S0 0, R0, L0, AP0)from a timed automaton TA = (S, S0, R, L, AP, C, λ, τ) 11.4. CONCLUSIONS 111 where: S∗=S0/≡is a finite set of states obtained by the equivalence classes of Eval(C)and S, S∗ 0={[s0]|s0∈S0 0}is the set of initial states, R∗⊆S∗×S∗is a total transition relation between states defined as: •Adiscrete transition (hsi, αi,hsj, α0i)such that αis an equivalence class and ri= (si, sj)∈R. The equivalence class αis not a boundary class. The equivalence class α0is the equivalence class αwith some clock that are reset. The clocks λ(ri)are reset, i.e., ν(x)=0,∀x∈λ(ri)and ∀ν∈α. The clock valuation νsatisfies a firing condition τ(ri). •Adelay transition (hsi, αi,hsi, succ(α)i)where αis not an end class. L∗(hs, [ν]i) = L(s)∪ {g∈AP∗\AP |[ν]g}is the labeling function. AP∗=AP ∪CC(C)∪CC(φ)with CC(C)the clock constraints of the timed automaton TA and CC(φ)the clock constraints introduced by the formula φ. Therefore, the verification of a TCTL formula φin a timed automaton TA is reduced to the verification of a CTL formula φ0in RTS(TA,φ). Once the region transition system is obtained, a CTL-like model checking algorithm is applied. Thus, the complexity of verifying a TCTL formula φagainst a timed automaton TA is linear in |φ|and (|S|+|R|), with |φ|the number of logical connectives and temporal operators of the formula, and |S|and |R|the number of states and transitions in the region transition system RTS respectively. More generally, the complexity is Θ((|S|+|R|)∗ |φ|). The commercial software UPPAAL is an example of a model checking tool that implements these algorithms and it is prepared for managing timed transition systems and logics. The open source NuSMV and PRISM also include variants of these logics and systems. In fact, NuSMV assumes that each transition of the transition system consumes a clock tick. Finally, PRISM implicitly embeds the notion of time in probabilistic timed automata and Markov chains. 11.4. Conclusions The main contribution of this chapter has been the study of phylogenetic properties that need the use of explicit time in the specifications. Firstly, we have motivated the necessity of this extension with real examples; for instance, the inspection of the date of appearance of certain phenotypes in the branches of the 112 CHAPTER 11. TIMED TRANSITION SYSTEMS AND LOGICS phylogenetic tree such as the peculiar traits of the Tibetan populations or woolly mammoths. Next, we have collected the requisites for augmenting the model checking framework and managing this new kind of restrictions in the specifications. To this end, we have introduced TCTL for solving these problems. TCTL is an existing temporal logic with real time capabilities: it provides a syntax, semantics, model checking algorithms and tools that fit with our requirements. After the presentation of the logic, we have shown how to translate a phylogenetic property to this particular notation. This approach shows the conceptual feasibility of our proposal for defining hypothesis with temporal restrictions and distances. Finally, we have presented the algorithms and software tools that can process this kind of formulas. The estimated theoretical complexity of verifying this type of specifications is Θ((|S|+|R|)∗ |φ|), with |φ|the number of logical connectives and temporal operators of the formula, and |S|and |R|the number of states and transitions in the region transition system RTS associated to the phylogenetic tree respectively. The experimentation with this framework is delayed to further chapters. In particular, Chapter 14 and Chapter 15 evaluate phylogenetic properties that mix probabilities with temporal constraints, which poses a more complex, versatile and realistic framework for the analysis of phylogenies. The time and memory consumption observed in the experiments over there places an upper bound for the properties studied here. Following to this chapter, Chapter 12 continues with the analysis of temporal distances between sets of tree states but from a complementary perspective: now, the output result of a phylogenetic property is an integer indicating a distance. Chapter 12 Computing Distances Between Symbolic Objects Time is the longest distance between two places. —Tennessee Williams, The Glass Menagerie 12.1. Introduction The result of some phylogenetic properties is not necessarily constrained to a boolean result. As well as checking the validity of a formula with qualitative model checking, biologists often need quantitative information as output of an evaluation. For example, the increasing number of inferred phylogenies and procedures for computing them leads to the necessity of incorporating comparative methods that help to analyze the structural differences of the trees. These comparative methods are based on distances between sets of states, the number of movements for converting one tree in another, or more generic topological measures (tree balance and asymmetry) involving the calculation of numerical values. In fact, the use of topological measures allows the comparison of heuristics for building phylogenies as well as the discrimination of trees violating phylogenetic restrictions [130, 158, 159]. Thus, the calculation of quantitative information is a key question for comparing and discarding trees. The identification and selection of the elements over which we will evaluate an arithmetic function comprises the first step. Here, the aim of this chapter is twofold. First of all, we focus on the symbolic manipulation of phylogenetic objects: we try to characterize a particular set of states of the tree satisfying certain restrictions imposed by a formula of temporal logic. Later, we show that the model checking framework supports the manipulation of these objects and the computation of distances among them. Previously in 113 114 CHAPTER 12. DISTANCES BETWEEN SYMBOLIC OBJECTS this thesis, the analysis of phylogenies with model checking was centered primarily on the specification and evaluation of boolean hypothesis over the whole tree. Now, we focus on the characterization of a particular set of tree states satisfying certain restrictions imposed by a formula of temporal logic. For instance, the set of nodes of the phylogeny belonging to a clade represents one of these objects. The computation of the distance between two clades is important for discerning if the tree is coherent or aberrant with respect to a hypothesis of the phylogenetists. Obviously, the inclusion of numerical restrictions in the computations requires the incorporation of quantitative capabilities in current temporal logics and model checking algorithms. The presentation of the TCTL logic (Chapter 11) provides a syntax and semantics for the definition and inquiry of properties with explicit integer values in the specifications. But the distances introduced in the specifications can also be discovered by the model checking algorithms in some cases. Following the same classification introduced for the labeling of the branches with quantitative information, the numerical output of model checking can be grouped in distance (or temporal) values and probabilities. In this chapter, we focus exclusively on returning temporal distances between phylogenetic objects defined symbolically as result of the model checking process. We show how some common phylogenetic metrics can be formulated with temporal logics. Alternatively to all the measures involving time and distances, we center our attention for probabilities in Chapters 14–15. There, we initially study the truth of a boolean property over a phylogenetic tree with a particular instance of the probability threshold. Most of the usual model checking tools are also capable of automatically infer the upper (maximum) and lower (minimum) probability bounds. 12.2. Returning Time Distances as Output of the Model Checking Procedure Firsts extensions of classic model checking algorithms cater for the computation of quantitative bounds. Normally, they count the maximum or minimum number of (timed) steps between groups of elements, considered as sets of states characterized by a logical formula [35]. These modifications in the logic result in the addition of MAX and MIN operators to the repertory of instructions of common model checkers such as NuSMV, which returns numerical values as output of their evaluation. Those operators belong to the syntax and semantics of the RTCTL logic, one of the firsts time extensions to CTL that assumes a weight of one time step in each branch of the tree. As root, common ancestors, leaves and clades are easily characterized sym- 12.2. RETURNING TIME DISTANCES 115 bolically with a proper labeling of the states and querying logic (see Chapter 4, Section 4.2), this feature provides a precise way for the computation of branch lengths and distances. A clear example of application corresponds to the calculation of the chronological date of divergence between populations in the Tibetan Plateau, where each one is characterized by its particular polymorphisms. Moreover, most tree (im)balance and asymmetry metrics for testing the quality of a phylogeny use topological distances between the root and subtrees or leaves. Sometimes, they discriminate the path lengths of left and right descendants (whose direction is marked when labeling the states of the tree) [136, Table 2]. For instance, the computation of the mean topological distance Mfrom the root to the leaves is: M=1 NΣi∈LMi where N=|S|is the number of states in the tree (Kripke structure), Lis the set of leaves and Miis the topological distance, i.e., the number of intermediate nodes between the ith leaf and the root. In the particular case of acyclic directed graphs such as trees, a pair of states are connected with a single path and thus the MAX/MIN distance between nodes is identical. Any of the next two expressions in NuSMV will return the same value for the length between the ith leaf and the root: COMPUTE NAME Mi := MAX [id=root, id=i] COMPUTE NAME Mi := MIN [id=root, id=i] Storing the partial results of each Mi and applying the correspondent addition to all the elements, we obtain the aforementioned metric. The mean topological distance is used for the calculation of the statistical variance or in the next metric: Σi∈LMi/2Mi As the range of output values of Mi is well-determined by the number of nodes of the tree, the balance of the tree is known by the examination of this result. Additional metrics are defined in terms of maximum distances between the root and one of the leaves: Σi∈IZ−1 i where Iis the set of internal nodes. The Zivalues are computed through: COMPUTE NAME Zi := MAX [(id=root) & !terminal, terminal] 116 CHAPTER 12. DISTANCES BETWEEN SYMBOLIC OBJECTS The tree states are enriched with boolean variables in order to identify internal and terminal nodes. Hence, it is patent that distance-based phylogenetic inference methods can take advantage of this particularities and metrics for the construction of a tree [176]. The previous approach is motivated by systems modeled with timed automaton and discrete time logics where the branches are labeled with unit time stamps. Further extensions should focus on the computation of MAX and MIN paths in branches tagged with more flexible time stamps or lengths, i.e., continuous time with t≥1. In addition, new and more complex quantitative functions can be applied over the result set characterized by a symbolic function. Powerful quantitative model checking tools such as PRISM allow us to count the number of states satisfying a property or relation, the average of a numerical label in those states or the maximum/minimum value. Supplementary balance metrics may take advantage of these functions. For example, the metrics involving the count of leaves in left (lj) and right (rj) subtrees descendant of the ancestor j: 2 (N−1)(N−2)Σj∈I|rj−lj| PRISM can evaluate the equation using the subsequent steps. Figure 12.1 represents the phylogenetic tree depicted in Figure 2.1 expressed in the notation of the PRISM model checker. Each state is labeled with an identifier and two boolean values: the first one tells if it is a leave or internal node, and the second one tells if the node is the left or right descendant of his parent in case of a binary tree. The phylogenetic tree is defined backwards: from the leaves to the top, plus a self-loop in the root. The name of the root is changed to X1 in order to avoid collisions with the reserved word X(next operator). The counting of the terminal leaves belonging to the left subtree of the root is achieved executing the following expressions. The first step returns the identifier of the left son of the root, and the second step calculates the number of terminal leaves that has the previous node as its root; i.e, the leaves belonging to the left subtree. The computation of the number of leaves in the right subtree is symmetric. "direct_left_son": filter(print, left & P>=1 [X id=X1]) filter(count, P>=1 [F "direct_left_son"], terminal) Additionally to these metrics, the direct comparison of phylogenetic trees is a key question. Robinson-Foulds [146] is one of the most common and popular metrics in comparisons of (un)rooted trees, but there are many more [102]. Another phylogenetic comparative method called SPR, that calculates the distance of two trees as the minimal number of moves that transforms one tree into the other, is Probabilities play an important role for the study of phylogenetic properties. They are applied in two complementary domains. From the point of view of the topological coherence, the additional labels of the tree mark the statistical consistency of the branches. Probabilities tell the number of bootstrapped or inferred trees that support that specific shape. There exists a range of statistical methods that compute these support values, which are usually expressed with percentages between [0,1]. One of the most common processes for scoring and comparing trees numerically is based on the computation of maximum likelihoods estimations [69]. Not only do they play an important role for the comparison of trees through the evaluation of the tree coherence, but also probabilities are useful for showing the prevalence of genetic patterns inside populations. Omitting the concrete support values in the branches of the phylogeny, we can center entirely on the study of the labels in the tree states. Even assuming that all the branches are generated with equal probability, the likelihood of finding a particular label (mutation or phenotype) in a state is not homogeneously distributed in the tree and it will variate with respect to each section of the phylogeny. Hence, the calculation of percentages of satisfaction of a property inside a clade or along a path of the phylogeny enlarges the expressiveness during the verification. In any case, the necessity of extending the formal framework with probabilities for both circumstances is evident. Although they share common roots, the operations needed for the computation of the topological support slightly differ from the operations needed for the computation of probabilities in reachable states. The first problem usually evaluates continuous time equations for the elaboration of maximum likelihood estimations, while the second problem is simplified with the use of discrete time distances in the paths of the phylogeny. The application domain changes with the type of time considered for the computation of probabilities. That is the main reason why the content of probabilistic model checking is divided in two sections. The aim of Chapter 14 is to present the main concepts of probabilistic model checking with a probability function that depends on the discrete time associated to the branches and, later, increment the complexity of probabilistic model checking with continuous time in Chapter 15. In particular, the first chapter is devoted to analyze the variation of a state property in several zones of the tree and compute the probabilities of satisfaction, while the formalism introduced in the second chapter is applicable to the estimation of maximum likelihoods. 123 Chapter 14 Discrete Time Probabilistic Transition Systems and Logics The probability of life originating on Earth is no greater than the chance that a hurricane, sweeping through a scrapyard, would have the luck to assemble a Boeing 747. —Richard Dawkins, The God Delusion 14.1. Introduction Previously, we have seen that the phylogenies are occasionally enriched with time labels or weights in the edges. This knowledge is useful for learning complex properties about the evolution, for instance, the estimation of the temporal point of divergence between species [17] or the diaspora of human populations [36]. The extension of the phylogenetic properties in Table 4.1 with time and probabilities increases the expressivity of biological hypothesis. Take the following disease as a clarifying example. The lactose intolerance in adults is a chronic disease caused by the inhibition of the lactase gene after the breastfeeding and childhood. The inability for processing the milk and its derivations is not homogeneously distributed in the human population. While in some African pastoralist groups in North/East Africa and the northern cultures of Europe their stock breeding tradition and diet motivated an evolutionary adaptation to digest the milk (>70% of tolerance), the percentage of acceptation decreases in the rest of areas and ethnic groups [160, 92]. In addition, the phenotype in Europe and Africa appeared at a different epoch and the point mutations that regulate the activation of the lactase persistence are disparate [163, 97]. Some illustrative questions that we desire to ask to the phylogeny, and that are expressed below, 125 126 CHAPTER 14. DISCRETE TIME AND PROBABILITIES require the addition of time to the branches of a population tree. The time allows the estimation of the divergence points between individuals or mutations, while the probability of the lactose persistence in different zones is calculated through the study of the distribution of the point mutations that regulate the phenotype. The questions are: 1. What is the rate of lactase persistence in a population? i.e., do their members define a characteristic haplogroup? and in that case, 2. Which polymorphism, among the multiple activators and inhibitors of the lactase gene, is the most frequent over there? and finally, 3. When did this phenotype approximately start to be predominant? i.e., does this date mark an important event in the diet, culture or migration of that population? These questions ask about the time (dates) and probabilities (frequencies/rates) stored in the branches of the tree. Besides, the deductive process that answers the queries also needs the manipulation of quantitative information. Thus, we must introduce a logic, a transition system and a model checking algorithm capable of expressing and managing these kind of questions. The notion of time introduced here matches with the concept of evolutionary or chronological clock. This chapter is organized as follows. After this introduction, Section 14.2 presents the syntax of the discrete time probabilistic logic. Secondly, Section 14.3 introduces the model checking algorithm for PCTL. Next, Section 14.4 focuses on the existing tools that handle this kind of logic and data structures. Finally, Section 14.5 draws the conclusions of this chapter. 14.2. Discrete Time Probabilistic Logic and Structure In this section we are considering a phylogenetic tree enriched with numerical information that tells the probability of selecting a branch descending from an internal node. Therefore, we can answer questions like: what is the probability of reaching a set of states of the tree from the root? The logic and data structure defined here settle the basis for future updates and extensions for continuous time systems in Chapter 15. Stochastic systems generally use Markov chains as the underlying data structure that provides semantics to the verification process [106]. Discrete time Markov chains capture the essentials of probabilities between states of the tree and implicitly associates an unit time step to every transition of the system. 14.2. DISCRETE TIME PROBABILISTIC LOGIC AND STRUCTURE 127 Definition 14 (Discrete time Markov Chain).A discrete time Markov chain is a finite transition system represented by a tuple M= (S, S0,P, L), where: Sis a finite set of states, S0⊆Sis the set of initial states, P:S×S→[0,1] is the transition probability matrix that indicates the probability of outgoing from state sito a certain state sjsatisfying Σsj∈SP(si, sj)=1, and L:S→2AP is the labeling function that associates each state with the subset of atomic propositions that are true of it. A phylogenetic tree is assimilated to a discrete time Markov chain making the corresponding associations of states to the definition of phylogeny (Definition 3). The leaves are labeled with the genome information of the population or specie they represent, plus additional data when necessary. Each branch of the phylogeny is mapped to an element P(si, si+1)>0 of the transition probability matrix. This value gives the probability of making a transition from state sito state si+1 in one time step. By default, we assume that all the descendants sjof a state siare equiprobable (P(si, sj)=1/n with nthe number of successors), but this value can be adjusted. The terminal leaves modeled with self-loops in the Kripke structure of a branching-time phylogeny are represented in the transition probability matrix by a single transition going back to the same state with probability 1. For any set of infinite paths Π starting in the initial state s0, the subset Π(πn) selects the paths π∈Π whose prefix equals to the finite sequence πn=s0s1s2. . . sn of length n+ 1 states. The set of infinite sequences sharing the prefix πnhas probability Pr(Π(πn)) = PΠ(πn). The probability PΠ(πn) is calculated as the product of probabilities for each intermediate transition, except for paths with unitary length in which case n= 0, π0=s0and PΠ(π0) = PΠ(s0) = 1: PΠ(πn) = 1 if n = 0 P(s0, s1)·P(s1, s2)·. . . ·P(sn−1, sn) otherwise Bayesian model checking methods allow for analyzing stochastic systems [100]. Probabilistic CTL (PCTL) [150, 88, 13] and Continuous Stochastic Logic (CSLTA) [59] help to formulate conditions on a discrete or continuous time Markov chain, respectively. The properties are referred to state formulas (φ) or path formulas (Φ). Besides, they allow enriched queries of temporal formulas such as P∼λ(Φ). Given an initial state sand a comparison ∼∈ {<, ≤,=,≥, >}, the operator P∼λ(Φ) returns true if the probability for a set of paths satisfying Φ is ∼λ, with λ∈[0,1]. 128 CHAPTER 14. DISCRETE TIME AND PROBABILITIES Definition 15 (Probabilistic Computation Tree Logic).A temporal logic formula φis defined by the following grammar with the minimal set of operators, where p∈AP and k∈N∪ {∞}: φ::= true |p| ¬φ|φ∨φ|P∼λ[Φ] (14.1) Φ ::= Xφ|[φU≤kφ] The formulas are checked against a structure Mconsidering all infinite paths π∈Πfrom a certain state s0. Notice that M, s0φmeans that s0satisfies φ. The semantics of well-formed formulas is as follows (let π=s0s1s2. . .): M, s0p⇔p∈L(s0), M, s0¬φ⇔M, s02φ, M, s0φ∨ψ⇔M, s0φor M, s0ψ, M, s0P∼λ[Φ] ⇔Prob(M, s0,Φ) ∼λ, The calculation of the probability Prob(M, s0,Φ) requires the identification of the infinite paths πsatisfying the path formula M, π Φ: M, π Xφ⇔M, s1φ M, π [φU≤kψ]⇔ ∃0≤i≤k, ∀0≤j≤i: (M, siψ)∧(M, sjφ) This set, {π∈Π|M, π Φ}, can be obtained by the union of finitely many pairwise disjoint subsets Π(πn)by [106, Definition 3], each one characterized by the finite prefix πnof all infinite sequences of the set. Therefore, Prob(M, s0,Φ) = Pr{π∈Π|M, π Φ}= ΣπnPr(Π(πn)) computes the probability as the summation of probabilities in all possible prefixes πnby [106, Theorem 1]. These logics usually support timed transitions in the Uoperator. The notion of time in a Markov chain falls within the concept of state distances. Each state transition of the discrete time Markov chain involves an unit time step. A mapping between the chronological time and state distances allows the inference of the evolutionary speed in the branches of the phylogenetic tree. The computation of time and probabilities are embedded in the model checking algorithm. Timed variants of the modal operators Fand Gare obtained via Uas F∼cφ=trueU∼cφ and G∼cφ=¬F∼c¬φ. Instead of writing intervals explicitly, sometimes they are abbreviated with inequalities. For example, P≤0.5[Φ] denotes P[0,0.5] [Φ]. Far beyond the use of discrete or continuous time, the main difference between PCTL and CSLTA syntax is the substitution of the long-run operator Lin PCTL 14.3. ALGORITHM FOR PCTL MODEL CHECKING 129 by a steady-state operator Sof CSLTA. These advanced operators are unnecessary for the phylogenetic problems we are trying to solve, but they could be included in the future if necessary [106, 13]. By now, we can translate the questions presented in the motivation example of lactose into the PCTL syntax. In a phylogenetic tree, the tips correspond to individuals of disjoint populations whose states are tagged with their DNA and a boolean indicating if they are lactose (in)tolerant. The internal nodes of the inferred ancestors are labeled with their estimated DNA sequence and lactose phenotype as well. The following equation asks if there exists an ancestor (P>0) at distance 3 or above from the initial state (F≥3) that is the root of a population with lactase persistence over 70% (P≥0.7[F≥0lactose tolerant]). The members of a population, including the leaves and internal nodes, are reached through F≥0. P>0[F≥3(P≥0.7[F≥0lactose tolerant] )] (14.2) The outer restriction P>0[F≥3] corresponds to the question 3 of the motivation. It searches for an internal node from which the phenotype starts to be predominant after a certain date since the phylogenetic root. The inner formula P≥0.7[F≥0lactose tolerant] answers the question 1 about the rate of lactase persistence in a population. Finally, the addition of a genetic marker in this place inside the P≥0.7equation helps to investigate the relation between a polymorphism and phenotype (question 2). The evaluation of the formulas needs the algorithm introduced in the next section. 14.3. Algorithm for PCTL Model Checking The evaluation of formulas written in PCTL over discrete time Markov chains drastically differs from the evaluation of formulas written in CSLTA over continuous time Markov chains due to the divergent interpretation of time. This point clearly arises during the verification process, which adds semantics to the specified formulas. The model checking algorithms for managing and solving PCTL or CSLTA formulas in stochastic systems are mainly identical to those of classic model checking except for the resolution of P∼λ[φ], i.e., the next and until operators with probability thresholds. In short, the recursive algorithm of model checking incorporates the new sentence: Sat(P∼λ[Φ]) = {s∈S|Prob(M, s, Φ) ∼λ} P∼λ[Xφ]formula. In PCTL, the probability of satisfying the next operator requires the probabilities of the immediate transitions from s. It is resolved by: Prob(M, s, Xφ)=Σs0∈Sat(φ)P(s, s0) 130 CHAPTER 14. DISCRETE TIME AND PROBABILITIES P∼λ[ψU≤kφ]formula. The computation of the probability for the until operator depends on the value of k. For the case of k∈N, then Prob(M, s, ψU≤kφ) is equal to:    1 if s∈Sat(ψ) 0 if k= 0 or s∈Sat(¬φ∧ ¬ψ) Σs0∈SP(s, s0)·Prob(M, s0, ψU≤k−1φ) otherwise When k=∞, the until operator is analogous to the original until operator of CTL with semantics of infinite paths. That is, Prob(M, s, ψU≤∞φ) can be rewritten as Prob(M, s, ψUφ) and it equals to:    1 if s∈Sat(ψ) 0 if k= 0 or s∈Sat(¬φ∧ ¬ψ) Σs0∈SP(s, s0)·Prob(M, s0, ψUφ) otherwise The time complexity of verifying a PCTL formula φagainst a discrete time Markov chain is linear in |φ|and polynomial in size of S, with |φ|the number of logical connectives and temporal operators of the formula. More generally, the complexity is Θ(poly(size(S)) ∗tmax ∗ |φ|) where tmax is the maximal step bound of a path subformula ψ1U≤tψ2of φ, with tmax = 1 if it doesn’t contain any Usubformula. The evaluation of CSLTA formulas demands the upgrade of Markov chains structures for providing continuous time semantics. The computation of probabilities is assisted by a partial transformation of the continuous time problem to a discrete time problem. Hence, the operations realized by the model checking algorithms for the calculation of probabilities in PCTL are the basis for future extensions. The time complexity of verifying a PCTL formula also poses a lower bound to the complexity of analyzing a CSLTA formula. The introduction of these semantics and complexity are presented in Chapter 15. 14.4. Model Checking Tools and Experimentation PRISM [124] is a generic model checking tool capable of handling probabilistic and timed specifications over Markov chains. Among its basic conceptions, PRISM checks if the probability of reaching a set of satisfiable states is up or below a 14.4. MODEL CHECKING TOOLS AND EXPERIMENTATION 131 predefined threshold. There exists a considerable diversity of model checking tools with different performances and qualities [98]. Although the real performance depends on the particular structure of the model and specifications, PRISM offers Java portability, a powerful syntax for handling time and probabilities in models and specifications, and a good scientific community support. Besides, it is open source, which allows the modification and optimization of its code. The model checking tool requires two input files for the verification process: a first file with the description of the model, and a second file with the specification of the properties. A description of the phylogenetic tree (Kripke structure with atomic propositions) in PRISM syntax musts be provided by the user as first input for the model checker. To this end, we precompute a sequence alignment and a phylogenetic tree. The data set used for this experimentation is synthetic. With this data set, we try to cover the spectrum of small phylogenies and analyze the cost of the evaluation of the lactose property over there. We have created random phylogenetic trees of up to 1000 tips using a Yule backward model [156]. For each tree size (number of tips), we have generated ten random trees and calculated the harmonic mean time. The DNA sequences have 50 bases with an homogeneous distribution of nucleotides. The PRISM codification of the phylogeny follows the same idea presented for NuSMV in Section 5.3. Figure 14.1 shows the implementation of the branchingtime phylogenetic tree of Figure 2.1 in PRISM code. The main module describes the topology of the evolutionary tree, where the names of the tree nodes (taxa) are defined numerically (1,...,5). PRISM distinguishes between the current state and the next state using the quotation mark (’). The second part of the description consists of a function returning the DNA string associated to each node, which evolves in synchronization with the tree skeleton ([id1] tags in PRISM). The translation of the phylogenetic tree to the PRISM syntax has been performed automatically by a BioPerl script [153]. The script can be upgraded in order to include extra features such as the generation of multiple instances of phylogenetic trees, bootstrapping and so on. We have evaluated the lactose formulas introduced in the motivation but enriched for the detection of polymorphisms in the DNA. The underlying objective consists of the identification of a correlated evolution between the lactose tolerance and patterns in the genome. Other studies such as [92] use cultural information for discovering this coevolution and the influence of a milk-based diet. The utilization of phylogenetic comparative methods and regression techniques establishes the essentials of this approach. The probability threshold of the internal P≥x[F≥0seq[i] = j] ranges from x∈ [0.1,0.9], with i∈[1,50] the position where we search for the polymorphism and ja certain nucleotide. Figure 14.2 plots the time required for the computation of 132 CHAPTER 14. DISCRETE TIME AND PROBABILITIES 50 ×9 formulas corresponding to the expansion of iand xfor all the columns of the alignment and probability bounds. All tests have been run on a Intel Core 2 Duo E6750 @ 2.66 GHz with 8 GB RAM and Linux. PRISM performs well for the verification of the lactose formulas in small phylogenies because it follows a polynomial trend in time with respect to the number of tips. However, it requires the integration of new technologies and solutions to scale for larger phylogenetic trees and specifications. In fact, we desire to find the value of x,iand jparameters for which the verification of the equation returns true. The definition of patterns is a common procedure, which intuitively leads to parametric model checking [32]. Nonetheless, mining for knowledge without prior information requires a more or less thorough exploration of the structure, which can be combinatorial in some or all of its dimensions. Although inherently parallel, the exhaustive inspection of potential solutions involves the test of large sets of formulas and an intensive use of the topology and information of each state. The application of parametric model checking for model exploration is a future extension that will increase the potential of our framework. Additionally, the phenotype is regulated by several polymorphisms and the likelihood of finding a certain nucleotide depends on the selected path. For those reasons, the model checking process demands the introduction of operators that automatically compute the intolerance rate and the probability of reaching a certain nucleotide. The calculation of parameters representing probabilistic bounds in the specifications are already managed by current model checking tools. The computational capabilities of PRISM grants the inference of maximum (minimum) probabilities of satisfying a property. 14.5. Conclusions In this chapter we have motivated the extension for the analysis of phylogenies via model checking using quantitative information. We have enriched the phylogenetic tree with respect to the tree in Chapter 11, where we only use explicit time. Here, we have also proposed the inclusion of probabilities in the branches of the tree because of its natural interpretation in the phylogeny as the main novelty. In particular, we have presented a phylogenetic example based on the lactose (in)tolerance that needs these kind of quantitative information. To this end, we have introduced an extended logic and data structure adapted for probabilities and time together with the algorithms and computations for managing them. Our first goal has been the increase of the logical capabilities for querying about the date of appearance and degree of distribution of mutations and phenotypes. Next, we have experimented with synthetic data in order to prove the feasibility of our approach with existing probabilistic model checking tools. PRISM is a