Graphics Processing Unit–Enhanced Genetic Algorithms for Solving the Temporal Dynamics of Gene Regulatory Networks

https://doi.org/10.1177/1176934318767889 Creative Commons Non Commercial CC BY-NC: This article is distributed under the terms of the Creative Commons Attribution-NonCommercial 4.0 License (http://www.creativecommons.org/licenses/by-nc/4.0/) which permits non-commercial use, reproduction and distribution of the work without further permission provided the original work is attributed as specified on the SAGE and Open Access pages (https://us.sagepub.com/en-us/nam/open-access-at-sage). Evolutionary Bioinformatics Volume 14: 1–16 © The Author(s) 2018 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/1176934318767889 Introduction The understanding of gene regulatory networks (GRNs)1,2 is one of the main current challenges in biology. A crucial problem is the determination of laws that govern the temporal dynamics of a biologic system from a particular state to a later one, through activating and repressing actions of some genes or their products (proteins) over other ones in transcription and translation processes. The identification of these actions is very important for developmental biology, for determination of the causes of certain diseases, and for the signaling of therapeutic targets. The set of genes and their interactions, which determine a particular behavior, can be represented by a directed network in which the nodes are genes/proteins and the links represent the activating and repressing actions of some genes over other ones. Based on experimental data from 2 different network states, the dynamical processes that transform from one state into another through evolution rules can be inferred and their characteristic parameters can be evaluated.3,4 As experimental data are often unclear or incomplete, this information can also help experimentalists in the search for new interactions or in the discard of false positives. One of the most used inference method consists in using Boolean rules, characterizing node states by 2 binary values (activated/not activated) that can be marked as 1 and 0, and links associated to actions by activator/repressor do not exist (which can be marked as +1, –1, and 0). A set of possible simple rules can be established, which determines the state of a node at any given time by the nature of the incoming links to that node at a previous time. The state of the network nodes can be refreshed synchronously or asynchronously. Given an initial and a final state, the objective is to determine which rules and interactions are applicable to each node, with a particular configuration of links, to reach the final state (or a state very close to it) in a number of time steps, starting from the initial state. In other words, we would like to find the network links and the evolution rules of the network nodes. In the case that the result is not unique, the results can be interpreted with biologic criteria or by experiments that confirm one option against the other Graphics Processing Unit–Enhanced Genetic Algorithms for Solving the Temporal Dynamics of Gene Regulatory Networks Raúl García-Calvo1, JL Guisado1, Fernando Diaz-del-Rio1, Antonio Córdoba2 and Francisco Jiménez-Morales2 1Department of Computer Architecture and Technology, University of Seville, Seville, Spain. 2Department of Condensed Matter Physics, University of Seville, Seville, Spain. ABSTRACT: Understanding the regulation of gene expression is one of the key problems in current biology. A promising method for that purpose is the determination of the temporal dynamics between known initial and ending network states, by using simple acting rules. The huge amount of rule combinations and the nonlinear inherent nature of the problem make genetic algorithms an excellent candidate for finding optimal solutions. As this is a computationally intensive problem that needs long runtimes in conventional architectures for realistic network sizes, it is fundamental to accelerate this task. In this article, we study how to develop efficient parallel implementations of this method for the fine-grained parallel architecture of graphics processing units (GPUs) using the compute unified device architecture (CUDA) platform. An exhaustive and methodical study of various parallel genetic algorithm schemes—master-slave, island, cellular, and hybrid models, and various individual selection methods (roulette, elitist)—is carried out for this problem. Several procedures that optimize the use of the GPU’s resources are presented. We conclude that the implementation that produces better results (both from the performance and the genetic algorithm fitness perspectives) is simulating a few thousands of individuals grouped in a few islands using elitist selection. This model comprises 2 mighty factors for discovering the best solutions: finding good individuals in a short number of generations, and introducing genetic diversity via a relatively frequent and numerous migration. As a result, we have even found the optimal solution for the analyzed gene regulatory network (GRN). In addition, a comparative study of the performance obtained by the different parallel implementations on GPU versus a sequential application on CPU is carried out. In our tests, a multifold speedup was obtained for our optimized parallel implementation of the method on medium class GPU over an equivalent sequential single-core implementation running on a recent Intel i7 CPU. This work can provide useful guidance to researchers in biology, medicine, or bioinformatics in how to take advantage of the parallelization on massively parallel devices and GPUs to apply novel metaheuristic algorithms powered by nature for real-world applications (like the method to solve the temporal dynamics of GRNs). KEywoRDS: Gene regulatory networks, evolutionary computing, parallel genetic algorithms, GPU RECEIVED: November 23, 2017. ACCEPTED: February 28, 2018. TyPE: Original Research FUNDING: The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the following research projects of the Ministerio de Economía, Industria y Competitividad (MINECO), and the “Agencia Estatal de Investigación (AEI)” of Spain, co-financed by FEDER funds (EU): MABICAP (Bio-inspired machines on High Performance Computing platforms:a multidisciplinary approach, TIN2017-89842P), COFNET (Event-based Cognitive Visual and Auditory Sensory Fusion, TEC2016-77785-P) and TOP4COG (Topological Recognition of 4D Digital Images via HSF model, MTM2016-81030-P (AEI/ FEDER,UE)). DEClARATIoN oF CoNFlICTING INTERESTS: The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. CoRRESPoNDING AUTHoR: José Luis Guisado, Department of Computer Architecture and Technology, University of Seville, E.T.S. Ingeniería Informática, Avda. Reina Mercedes s/n, 41012 Sevilla, Spain. Email: [email protected] 767889EVB0010.1177/1176934318767889Evolutionary BioinformaticsGarcía-Calvo et al review-article2018 2 Evolutionary Bioinformatics ones. To find the structure of the network and its evolution rules, a heuristic optimization method is applied. A fitness function is calculated by comparing the reached state with the preestablished final state. The most frequently applied methods are evolutionary algorithms,5,6 a family of optimization methods that are inspired by the principles of natural evolution to find solutions to hard optimization problems. The method introduced in Aguilar-Hidalgo etal3 can infer the dynamics of GRNs between known initial and ending network states by simple acting rules, using a genetic algorithm (GA). The method was also reviewed and its applicability to real biologic problems was further discussed in Aguilar-Hidalgo et al.4 However, in those works, only a simply sequential version of the algorithm was used, thus needing several hours for a small testing GRN. Even for our optimized sequential implementation using a compiled language like C++ with all the optimizations enabled (used for comparison with a parallel version in this article), a runtime of the order of 30 min is needed, for example, to run a single-threaded version of the algorithm for a GRN including 32 genes with a GA population of 30 000 individuals on a high-end Intel i7 CPU. Moreover, this runtime increases rapidly with the number of genes of the GRN. Therefore, it is very convenient to design and implement efficient parallel algorithms for this task. Genetic algorithms have an inherently parallel nature that makes them well suited to implementation in massively parallel hardware like graphics processing units (GPUs), but they have several sequential parts that are difficult to deal with, like the crossover or the selection procedure. Thus, it is important to fine-tune the parallel implementation to benefit from all the potential computational power of the GPU. To allow the practical application of this method to real-life problems (often involving large GRNs), we present in this work a comprehensive study on how to develop efficient parallel implementations of this method for the fine-grained parallel architecture of GPUs. To this end, we evaluate both performance and effectiveness for a wide range of possible elections for parallel GA schemes (master-slave, island, cellular, and hybrid models) and selection methods (roulette, elitist), and we suggest procedures to optimize the use of the GPU’s resources. Main contributions Because of the intricacy of the GRN dynamics, we have proceeded in a systematic manner to unveil what the most promising GA and the most efficient GPU configurations may be. All along the article, our conclusions are described and justified with experimental data. As a summary, which can be a guide for practitioners and experimentalists to solve the temporal dynamics of GRNs, the working flow of our experiments was as follows: 1. Studying the transformation rules that can be applied to each network state and finding a universal encoding for them. 2. Defining a proper fitness function based on the number of differences in network node and network link values, between the current population individual and the target solution. 3. Implementing and testing most of the parallel versions of GAs, bearing in mind 2 main groups: nonstructured master-slave models and parallel structured GA models (island, cellular, and hybrid). 4. Confronting the behavior of 2 dissimilar selection methods: roulette and elitist. 5. Studying in depth the best model found. In particular, we have discovered a large variation in the results for the island model when varying the migration frequency (number of generations between each migration) and the migration rate (number of individuals migrated). 6. Determining whether the found fitness is near to the optimal or not. In fact, we have found the minimum possible fitness along our experiments. Although there is much general advice when programming for GPGPUs (General-Purpose Computing on Graphics Processing Units), for the specific case of achieving maximum performance in the study of GA for GRN some important aspects have been unveiled all along this work. Several key points to achieve maximum performance have been found accordingly. In particular, we have determined a different thread pairing for each genetic operator as a parallelization strategy that GAs may benefit from. The complete list of those key points is as follows: 1. First, the proper selection of the number, type, and size of the data objects for the GPGPU is of essential priority. 2. Determination of the best thread pairing for each genetic operator (or any other computation). The most extended and intuitive thread selection in the current literature (pairing individuals with threads) does not usually produce the better results. 3. Reorganization of memory accesses to favor coalescence in the memory access patterns. 4. Careful estimate of GPU shared memory consumption so as the maximum number of threads could fit in it. This estimation usually determines the maximum number of individuals for some models of the GA, such as the cellular model. 5. Management of any type of bifurcation (branch, loops, etc.). On the whole, we have found that it is especially important to avoid thread divergence when applying the GRN rules. The execution time of the algorithm has been reduced by applying development techniques which optimize the use of resources of the GPU with compute unified device archit (CUDA). A comparative study of the performance obtained by the different versions of the parallel implementation on GPU against a sequential implementation on a CPU has been made. García-Calvo et al 3 The proposed parallel algorithm has been implemented into the software cuEVOGENET and is available online at https:// github.com/rxp90/cuEVOGENET. Related Work Evolutionary algorithms are based on creating a population of possible solutions of a problem and applying procedures that mimic natural evolution to find good candidate solutions within that population: creation of diversity, mechanism of selection, and genetic inheritance. Thus, evolutionary algorithms explore the search space to locate good solutions of the problem. Genetic algorithms7,8 are a class of evolutionary algorithms in which a candidate solution to a given problem is encoded as a bit string, called individual or chromosome. A population of such strings is chosen at random and evolves over a number of generations under procedures inspired in natural evolution to find good candidate solutions: selection according to a fitness function, crossover to produce new offspring, and random mutation to increase diversity. Genetic algorithms and other types of evolutionary algorithms have been successfully applied to various bioinformatics problems, such as inferring GRNs from gene expression profile data,9,10 identifying motifs in DNA sequences,11,12 or generating DNA strands for DNA computing.13 Many of the GA operations can be executed in parallel in different processors or processing elements with a high computation to communication ratio. Therefore, GAs are good candidates to be executed in parallel and distributed computers with good speedup respect to sequential executions. This parallelization is very often necessary to obtain solutions of good quality in a reasonable time. Even though GAs are based on simple mechanisms, parallel GAs are complex nonlinear algorithms controlled by many parameters, so that it is difficult to preview which are the best choices for a good efficiency. Different models of parallelization of GA have been studied. A good early overview of techniques employed to implement GAs in parallel is the work by Alba and Troya.14 Cantú-Paz15 also presented a classification of parallel GAs and made a systematic study of how the different parameter values affect to the quality of solutions and to the efficiency of the parallel implementation. A recent survey of the state of the art in distributed evolutionary algorithms (including and up-to-date references) is offered by Gong etal.16 Moreover, this work analyzes in detail the different models of parallel evolutionary algorithms and their characteristics. In addition to the runtime reduction, another advantage of the parallel implementation of evolutionary algorithms is that it can add intrinsic fault-tolerance mechanisms to the optimization algorithm, as shown by González etal17 for the master-slave model and by Hidalgo etal18 for the island model. Parallel GAs have been implemented on GPU for many test problems using various parallelization models, demonstrating that good results can be obtained specially for large populations. Some examples are Luo and Liu,19 who employed a cellular model; Li etal,20 who used a master-slave model; and Luong etal,21 where island model was preferred. Some GPU implementations of GAs are also reported for real-life biologic problems. For instance, Ben-Shalom etal22 used a master-slave model for a computational neuroscience problem consisting in tuning the parameters of kinetics models for voltage-dependent ion channels of neurons, obtaining a considerable speedup for that problem as compared with a small Linux cluster. At last, a good review paper about parallel implementation of evolutionary algorithms is the work by Arenas etal,23 which introduces the questions posed by the implementation of this kind of algorithms on GPU and presents a survey of solutions found in the literature. In recent years, increasing attention has been paid to evolving GRNs. In consonance, a special issue of the journal BioSystems was devoted to this topic in 2009,24 identifying the 3 areas that have been more studied since then. The first area is synthesis of GRNs using evolutionary algorithms. Esmaeili and Jacob25 employed evolutionary algorithms to discover regulatory networks that optimize multiple stability indicators, including network sensitivity, cyclic length of the attractors, and number of attractors. Nicolau and Schoenauer26 used evolution to find regulatory networks with specific patterns of connectivity, such as scale-free networks. In the same line and most recently, Noman etal27 employed an evolutionary algorithm to find GRNs that include robust network modules, using a Monte Carlo method to quantify topologic robustness. They verified their algorithm by identifying GRNs showing 2 classic behaviors: oscillation and bistability. The second area is not related to this work (the reconstruction of biologic genetic networks using computational gene regulatory models and experimental data), so these references are not included here. Finally, the third area is the application of gene regulatory models and artificial development to solving engineering problems and to build artificial simulated organisms intended to study the working principles of natural systems. In the last decades, artificial GRNs, which are dynamical systems inspired in their natural counterparts, have been employed with those purposes. Evolutionary algorithms play a crucial role in finding good solutions of these artificial GRNs. Eggenberger28 proposed an artificial evolutionary system to evolve 3-dimensional shapes of simulated multicellular organisms. Banzhaf29 introduced an artificial regulatory network and showed that it can reproduce phenomena found in natural ones, such as heterochrony (a variation in timing of expression) or different dynamical behaviors on protein concentration, ruled by changes in the genome. Knabe etal30 proposed an artificial GRN to model biologic clocks and showed that it captures many of their characteristic properties. Guo etal31 proposed a GRN-based algorithm to make an arbitrary number of robots to self-organize into different predefined shapes and even to self-reorganize adaptively under dynamic environments without the need of a centralized control. They employed a multiobjective GA to tune the parameters in the regulatory model to simultaneously minimize 2 objectives: the 4 Evolutionary Bioinformatics total traveling distance of all the agents and the system convergence time. Joachimczak and Wróbel32 used artificial GRNs evolved by a GA to generate or process signals in which information is encoded in chemical pulses, showing their capability to perform tasks such as doubling the frequency or pulse length of the input signal, acting as a low pass frequency filter, or doubling the number of input pulses. More recently, Cussat-Blanc etal33 employed a special type of GA (with augmenting topologies) to evolve artificial GRNs that were used to control agents performing different problem domain tasks. GA for the Temporal Dynamics of GRNs In this section, we present a description of the GA that we have employed. A pure sequential, elementary implementation of it was introduced by Aguilar-Hidalgo etal.3 We refer the reader to this reference for more details on the basics of the algorithm. Although this previous work suggested a modelization for this kind of GRN, it did not give an insight in the most efficient GA model nor in reaching the best solution for this network. Maybe this was due to its elevated computing times (several hours for a trial), which made prohibitive the testing of the different possibilities: GA models, parameter adjusting, best number of individuals, mutation probabilities, and so on. Thanks to the chromosome encoding and the GPU-oriented parallelization methods employed here, we verify that the proposed GRN temporal dynamics model provides complete efficiency for the solution of this kind of problem. Indeed, we investigate all of these possibilities and find the implementation that produces better results (both from the performance and the GA fitness perspectives). An important result of our exhaustive study is that a moderate population size is enough to get to very good fitness values. This allows us to compute many of the GA possibilities in less than a minute and finally to find a minimal fitness value. Indeed, we have demonstrated that this value corresponds to the best solution for the analyzed GRN, impossible to diminish. Gene regulatory networks A GRN is modeled, using a Boolean network model,34 as a directed graph GVE=( ,) , where V is the set of vertices or nodes, which represent genes or proteins, and E is the set of edges on V , which represent regulatory interactions between genes/proteins. An edge between Vi and Vj symbolizes a dependence between them. The state of a node can be described by a Boolean variable expressing that it is active (1) or inactive (0), and thus that its products are present or absent. The graph will be annotated so as to reflect the nature of the dependencies (eg, activator, repressor). GRN temporal dynamics model The time evolution of a GRN starts from an initial state A and reaches a final state B . It is expected that the transformation dynamics that turns network state A into state B will be ruled by a certain rule set. The objective of the GA is to find a rule set (1 rule for each network node) that, applied to an initial state ′ A which is close to A and differs from it on a number of links, transforms it into a final state ( B *) that is close to the desired final state B , as shown in Figure 1. Each chromosome represents a candidate solution of our problem, including the set of links of state B * and a set of rules. ′ A is the state formed by the same node configuration as A and the same set of rules as B *. The application of this set of rules to state ′ A (1 rule for each network node) transforms it into a network state corresponding to B *. A pseudocode of the GA is shown in Figure 2. At the start, a population of solutions is generated. Each individual is evaluated according to a fitness function. After the generation of the initial population, evolution is allowed to proceed in an iterative manner. In each iteration or generation, following the biologic evolution analogy, several genetic variation operators (crossover, mutation) are applied to introduce diversity in the population. Each resulting population individual is evaluated according to the fitness function. Then a selection procedure is applied, in which the population for the next generation is formed by a combination of some individuals of the previous generation and the offsprings resulting from the application of the variation operators. The procedure is iterated until some stop condition is fulfilled: a predefined number of generations is completed or the fitness value of some individual surpasses a predefined value. Initial population The initial population is composed of individuals whose links part is generated by making random permutations of the original network links, that is, by permuting a number of link elements randomly chosen. In this way, a new individual with Figure 1. Evolutive transformation dynamics of the GRN. The objective of the genetic algorithm is to find a set of transformation rules that brings from ′ A(which is close to the initial state A and differs from it only on a number of links) to a final state B* close to the known final state B. García-Calvo et al 5 different links is obtained, but the connectivity of the network number of links of the individual is maintained. The rules part of each individual is generated randomly. Solution encoding The GA chromosome, which represents a possible solution of the problem, is formed by the model parameters whose configuration is to be investigated: the links among the nodes and the rules applied to each node. The chromosome used in our study is composed of 32 nodes, but for clarity we have shown in Figure 3 the chromosome of an individual with only 3 nodes. The first part is composed of the links that each node has, and the second part includes the rules applied to each of the 3 nodes. The possible values of a link can be +1 (if the link is an activator), –1 (if the link is a repressor), or 0 (if there is no link). A numeric value represents the rule to be applied to each node (see next section). Transformation rules The possible rules that can be applied to each network state are as follows: 1. Rule of the most: the final state of the node is determined by the action of most links directed to it. If most of them are activators it will be activated; if most of them are repressors, it will be deactivated. If the number of activators is equal to the number of repressors, the state of the node will remain the same. 2. Absolute repressor: the presence of just 1 repressor link is enough to deactivate the node. 3. Joint action of 2 activators: the node will be activated only if, at least, 2 input links are activators. 4. Joint action of 2 repressors: the node will be deactivated only if, at least, 2 input links are repressors. The set of rules is the specification of which one of these rules is applied to each particular network state in the GRN. Only links coming from active nodes are considered. Solution evaluation Each individual of the population is evaluated, and a score value is associated to it according to a fitness function, based on the number of differences in the values of the network nodes and the network links between the population individual ( B ) and the target solution ( B *). Hence, the fitness function φ is defined as, φλ λλ =(, ) (, )1(, ) (, ),0 1 dBB dBB dLL dLL * * * * maxmax    +− ()    ≤≤ (1) where dBB (, ) * is the distance between the values of the nodes in network states B and B *, and dLL (, ) * is the distance between the values of the links for B and B *, considered both distances as the number of differences between the nodes (or links) values. The better the individual, the lower its numerical fitness value. Here, λ is a regulatory parameter that determines the relative importance of each term (nodes or links) in the global calculation of the evaluation. As in the article where the model was introduced,3 the following value for the λ parameter has been chosen: λ = 0.9 (2) The reason for using a high λ value is that the amount of link variability found in these kinds of biologic problems is usually much higher than that of nodes. Genetic operators Crossover. In each iteration of the GA, just the best 2 individuals from the previous generation are kept, that is, elitism is applied. The rest of the individuals of the population are generated by applying crossover to couples of the remaining individuals. To select the parent individuals of each crossover operation, the roulette method is applied: all the individuals are capable of being selected as parents with a probability that is directly proportional to the normalized fitness value of the individual. Thus, individuals with a higher fitness value have a higher probability of generating new individuals and keeping part of their genetic information in the population. As a result, the fitness value of the average individual should increase with the number of generations. However, the crossover operations also introduce some variability so that the capability of exploring the search space is increased. As the chromosome is formed by 2 parts, rules and links, it is necessary to determine 2 related crossover points for each crossover operation. The links crossover point is determined randomly and the rules crossover point is determined according to the former election following the expression: Figure 2. Pseudocode of the sequential version of the genetic algorithm. Figure 3. Chromosome example representing a particular possible solution of the problem. 6 Evolutionary Bioinformatics Kint l N =1     + (3) where K is the rules crossover point, l is the links crossover point, and N is the number of nodes in the network. In this way, we calculate how many nodes are involved (ie, with all of their inputs being crossed) in the crossover and select also the rules attached to them. A low connectivity is frequently observed in biologic networks.35 Thus, to keep the connectivity of the network within low values, each new individual generated by crossover has to pass a connectivity test before being accepted as part of the new population. Only the individuals with a number of inputs per node smaller or equal to 7 are accepted. Otherwise, the individual is rejected and a new individual is generated by crossover. Mutation. After a new individual has been generated by crossover, it undergoes a mutation process in a probabilistic manner. There is a mutation probability for each one of the parts of a chromosome: µ1, the mutation probability for the links portion of the chromosome; and µ1, the mutation probability for the rules portion of the chromosome. Mutation is individually applied to each element of the chromosome according to its probability value. In the case that mutation has to be applied for an element, a new value for it is randomly chosen. All the possible final values after mutation have the same weight in this random choice: 1/2 for links (changing to 1 of the 2 remaining possibilities) and 1/3 for rules (changing to 1 of the 3 remaining possibilities). Selection/replacement. The population of generation t+1 is formed by the best 2 individuals from generation t and the individuals obtained by applying the crossover and mutation operators to the remaining individuals from generation t . Thus, the total number of individuals in the population is unchanged. GPU Architecture and Parallel Programming During the last decade, several reasons have pushed up GPGPUs as an excellent platform to conduct heavy tasks with a high degree of parallelism. One of the reasons for their success is its inherent general-purpose manycore-based implementation, in contrast to other less flexible architectures like Application- Specific Integrated Circuit (ASICs) or Field-Programmable Gate Array (FPGAs). Moreover, several additional advantages have thrusted their success for an ample set of scientific applications. First, the very low cost per processing unit, which is supposed to persist, since data center and PC graphics markets subsidize GPUs. Second, the annual growth of performance and memory bandwidth is predicted to remain very high.36 For example, many existing GPUs, with around 3000 floating-point arithmetic cores, can achieve a performance above 5 Tera FLoating point OPerations per Second (TFLOPS) and around 300 GB/s within 1 chip. One of the consequences of this advantage with respect to its competitors is that many top leading supercomputers are GPU-based machines.37 Besides, with the advent of new architecture capabilities of GPGPUs, and the expansion of tools, libraries, and hardware abstraction mechanisms, GPU software development is not now a tricky and bizarre programming task like it was in its former stages.38 Debugging and profiling tools like those from CUDA39 are more and more reliable. Previous reasons have pushed the scientific community to incorporate GPUs in several disciplines, and it seems that this tendency will continue for the next decade. Compute unified device architecture is the NVIDIA programming platform that has been chosen to implement current work. A GPGPU39 includes an array of streaming multiprocessors (SMs). Nowadays, each SM consists of 8, 16, or 32 floatingpoint Scalar Processors (SPs); 1, 2, or 4 Special Function Units; 1 or 2 multithreaded instruction units; and several pieces of memory. Each memory type is intended for a different use; thus, if the programmer were conscious of memory management, it could extract any drop of performance. Each GPGPU has nowadays 1 or 2 warp schedulers that select at any cycle several groups of threads for execution on each SM in a round-robin fashion. A warp is simply a group of 32 hardware-managed threads. As the number and types of threads are enormous, a 5-dimension organization is supported by CUDA, with 2 important levels: a grid contains blocks that must not be very coupled while every block contains a relatively short number of threads (usually 512 or 1024), which can cooperate deeper (eg, by sharing data through on-chip shared memory and fast synchronization primitives). If a thread in a warp issues a costly operation (like an external memory access), then the warp scheduler switches to a new warp, to hide the latency of the other thread. To effectively use the GPGPU resources, each thread should operate on different scalar data, with a certain memory access pattern. Implementation of the GRN Temporal Dynamics on GPUs Parallel-oriented GA models We have implemented different parallel versions of the GA: a nonstructured master-slave model and 3 parallel structured GA models—island, cellular, and hybrid model (Figure 4). In a nonstructured model, the GA uses a single population, and the GA procedures or operators (selection, crossover, mutation) are applied on them as a whole. In particular, all individuals can undergo crossover with one another (panmixia). In contrast, in a structured model, the population is divided in groups so that each individual can only interact with individuals within its group. Master-slave model. In this model, the GA is composed of a single monolithic (panmictic) population. To use the huge parallelism of the GPU, different operations of the algorithm can be executed independently for each individual of the García-Calvo et al 7 population on different processing elements of the GPU (see section “Optimizations for GPUs”). It is worth to mention that the implementation of the elitist selection has been targeted to reduce the GPU memory bandwidth consumption. Given the fitness values of an initial population (see Figure 5), sorting results are collected over their indexes. These indexes are then sorted by fitness values, so that adjacent couples are selected. The fourth row of Figure 5 shows an example of this selection, whereas its last row indicates which individuals (represented by their fitness values) are actually to be mated according to this criterion. Parallel island model. In the island model, also called “coarsegrained” or “distributed” GA model,15,40 the population is partitioned into a set of subpopulations or islands. An isolated GA is executed in each island and a number of individuals are copied to another island (migrations) periodically, as shown in Figure 4. This strategy has 2 advantages. First, a performance improvement can be obtained on a parallel computer. If we can execute each island concurrently on a different processor, the computing time of the algorithm can potentially be decreased substantially. Second, the presence of parallel populations and their interaction introduce a new dynamics that improve the numerical behavior and execution time of the algorithm by obtaining a higher genetic differentiation and a better sampling of the search space. A pseudocode of our parallel version of the GA for GRN temporal dynamics using the island model is shown in Figure 6. Different choices or parameter values are involved in the design of a parallel island GA: size and number of islands, migration topology (topology that interconnects the islands, for example, ring, star, or fully connected), migration frequency (number of generations between each migration), migration rate (number of individuals to migrate), and migration policy (selection of those individuals to migrate and those to be replaced in the destination island). As will be discussed in the section “Experimental results,” the choice of parameter values can affect the efficiency of the parallel algorithm. The methodology to implement migrations in the parallel island model is described in Figure 7. Information on islands was stored in the GPU memory. An elitist selection was carried out in each island. Then, the crossover operator is applied. To save computing time, the population was not sorted again after crossover: it is assumed that the crossing of good individuals will give rise also to good individuals. After that, the best n elements of each island are migrated to the destination island (being all the islands arranged in a ring topology). Parallel cellular model. In a cellular model (Figure 4), the whole population is arranged in a grid with a certain topology, and each individual can only mate and compete with other individuals included in its neighborhood. The neighborhoods overlap, so that individuals with good fitness value can propagate through the entire population. A cellular model is thus fine-grained and spatially structured. Figure 4. The 3 different models of parallel structured genetic algorithm that have been studied. Figure 5. Elitist selection implementation for the master-slave model. From top to bottom: fitness values of an initial population, increasing indexes, fitness sorting results expressed over the indexes, adjacent couples to be mated, and individuals to be mated according to this criterion. Figure 6. Pseudocode of the parallel version of the genetic algorithm following the island model. 8 Evolutionary Bioinformatics We have implemented a cellular parallel model in which elitist selection is used: each element is crossed with its best fitness neighbor as shown in Figure 8. Hybrid model. A hybrid model is composed of islands that include a cellular model inside each one (Figure 4). Hybrid models have not yet been used as frequently as island or cellular ones, because they involve a higher level of complexity, both in programming and in choosing the right values of the parameters of its more complex topologic structure. However, they can be good candidates to extract the potential power of parallel architectures with a highly hierarchical memory, such as GPUs. Optimizations for GPUs On the whole, GAs, as most bioinspired algorithms, are parallel in nature, which might make them good candidates for massively parallel hardware like GPUs.23,41,42 However, there are several sequential parts that are difficult to deal with (like the selection procedure). Hints about GPU optimization of these parts are unveiled in this section. One of the most important aspects that reduces time computation is executing everything directly at the GPU, thus avoiding transfer times between CPU and GPU. With current GPU architecture capabilities and software tools (mainly the all-purpose cuRAND, cuBLAS, and Thrust), this is perfectly possible and this policy has been followed all along this work. Some authors have even published an ordered points checklist for maximizing GPUs’ efficiency.43 In this respect, outstanding soft computing methods and excellent researchers meet with a profound challenge when trying to program in highly parallel platforms.44 On the other side, it can be easily understood that tuning GAs to be highly parallel and efficient over GPU platforms leads to implementations that are difficult to implement into other architectures. Although previous studies serve to depict a first scene of the main optimizations and challenges that a GA designer must face to when programming for GPGPUs, for the specific case of GRN other important aspects have been unveiled in this work. A lot of preliminary tests were carried out to discover the impact on performance of several programming decisions. In some occasions, we have detected that the most proper choice depends even on the GA model. During this evaluation, monitoring internal GPU performance through the Compute Visual Profiler45 for NVIDIA technology was used. Through the use of this tool, we have found several key points to achieve maximum performance. First, the proper selection of the number, type, and size of the data objects for the GPGPU is essential. In our case, GRNs contain lots of state information that consists on binary variables in most occasions. Although compressing this information in binary format would have supposed a reduction of total amount of consumed memory, a set of initial test did advise us that this compression was not beneficial from the performance point of view. The implementation of bit management operations was the main cause of this. As a consequence, we comprised most of the individual information, but each state value was stored in a byte. Special caution has been taken with shared memory consumption. The total amount of shared memory for the different parts of the GA has been carefully measured so as the maximum number of threads could fit in the 48 kB of the tested GPU (see section “Experimental results”). The second important decision is the proper launching of thousands or millions of very light threads. The most intuitive thread selection is pairing individuals with threads (all cited references here except one follow this method). Nevertheless, as some authors have pointed out, other possibilities can produce better results. For example, in Jaros,46 the author used a per-warp assignment, because in this case (solving the knapsack problem) warp granularity did not limit the size of individual (using appropriate mapping). In our case, we have stretched out a little more the thread assignment, mainly with the aim of avoiding GDRAM accesses by fitting the used memory with that of shared memory. This implies a different thread pairing for each genetic operator or other computations: •• Crossover. Each block crosses only 2 individuals. A reduction process is applied to the 512 threads of a block to compute if the connectivity of successors preserves the correct number of links according to GRN hypothesis. •• Mutation. A thread-individual pairing is preferred because mutation rules do not need individual communication. •• Selection/replacement. The pairing for this step must be carefully decided for each genetic model. Our results clearly reveal some pairings for each model, which can be extended for other GAs. Briefly, we can state that (a) if roulette is chosen, each individual has no relation with the others, and only a GDRAM access to the random numbers is necessary; thus, a thread per individual selection pairing is preferred; and (b) for elitist selection, CUDA libraries (like Thrust47) propose a sorting by key with a reduced number of threads. Moving pointers instead of moving the complete structure of individuals reduces the precious memory bandwidth, which yields to important benefit on CPU, which is even larger for GPU times. Figure 7. Description of the migration methodology used in the parallel island model. García-Calvo et al 9 •• Rule application. A node is paired with a thread. As the number of nodes of the analyzed GRN is 32, a warp contains an individual. For each node (thread), only that useful information (connected links and their state) is uploaded to the shared memory, thus reducing memory transferences. •• Fitness computation. Same pairing as rule application because similar reasons apply. Previous thread assignment and data partitioning to fit onchip shared memory reduce considerably main memory accesses. However, if this assignment wants to be preserved to keep previous benefits, some consequences on the maximum number of individuals appear for some GA models. Section “Study of fitness of cellular parallel model” gives details for this circumstance. The next decisive factor on performance is the management of any type of bifurcation (branch, loops, etc.). On the whole, thread divergence must be avoided for most cases, meanwhile coalesce memory accesses should be reinforced to maximize GPU utilization. There are several conditional structures in GRN dynamics where CUDA compiler is not able to introduce predication. These were detected using the NVIDIA Visual Profiler45 and considered of high priority by this tool. The most evident example is the choice of one of the 4 rules for each node (see section “Transformation rules”), when the rules are iterated 5 times for each individual. Other example is the selection of a random value after a mutation. These have been rewritten using lookup tables and logic bitwise operations, which introduce an important speedup. Another issue that must receive attention is the pattern memory access. The CUDA manual39 provides detailed algorithms to identify types of coalesced/uncoalesced memory accesses. Some reorganization of memory accesses was introduced in our codes to favor coalescing. Crossing operation requires memory accesses that may have a variable stride if not enough attention is paid. In our case, instead of taking the most evident decision, which is that any thread copies 2 consecutive nodes, the code was rewritten so that a thread copies 1 node and other with a distance of 512. This allows a stride equal to 1 between threads. At any piece of code, some speedup can be achieved. A simple test with 1 kernel shows different speedups for a block when evaluating 1024 individuals (using elitist model). Avoiding thread divergence when applying rules can achieve a 1.44× speedup, an acceleration of 1.66× is obtained when vectorizing loads when accessing individuals (using casting operators to read 1 int4 elements at once instead of 16 independent bytes), and finally, the vectorization of node link accesses when computing fitness values decreases time by a factor of 1.26. It is worth to mention that previous optimizations can be used in other parts of a GA, thus obtaining additional speedups. A relevant case occurred for migration of individuals when the island model is used (see section “Study of fitness of parallel island model”). Timing may vary significantly depending on the pairing between threads and individual information to be transferred. A comparison for 512 individuals (with a rate migration of 16 individuals each 16 generations) between 2 cases can give us an idea of the importance of this issue. In a first test, a thread was in charge of transmitting 4 bytes (1 int). In a second test, transferring was vectorized, so each thread transfers 16 bytes (1 int4). Because the resultant bandwidth is 10 times higher for the vectorized case, the execution is speeded up 3.5 times. Experimental Results A thorough study of various parallel implementations for the temporal dynamics of GRNs was carried out with a twofold aim: corroborating that GPUs are an excellent choice due to the huge achieved speedup with respect to a sequential CPU and selecting the best GA model (from fitness values evolution, scalability, and performance points of view). To achieve this goal, we have analyzed different parallel implementations for the GA which is the base of the studied algorithm: masterslave, island, cellular, and hybrid models. Moreover, we tested for GPUs the behavior of 2 different individual selection methods: roulette and elitist. Figure 8. An example showing how the best fitness neighboring element is chosen for each element in a 9-individual cellular population. Top: from left to right, first we depict the fitness value of 9 individuals arranged as a cellular model. On the next 3 diagrams, the Moore neighbors (considering a torus topology) for the individuals situated at (1, 1), (0, 0), and (2, 2) are marked with shadows. Bottom: left, coordinates of the best neighbors for the 9 individuals; right, a possible result after the crossover operation. 16 Evolutionary Bioinformatics 4. Aguilar-Hidalgo D, Lemos M, Córdoba A. Evolutionary dynamics in gene networks and inference algorithms. Computation. 2015;3:99–113. 5. Bäck T, Fogel DB, Michalewicz Z. Handbook of Evolutionary Computation. Bristol, UK: IOP Publishing; 1997. 6. Simon D. Evolutionary Optimization Algorithms. Hoboken, NJ: Wiley; 2013. 7. Goldberg DE. Genetic Algorithms in Search, Optimization, and Machine Learning. Reading, MA: Addison-Wesley; 1989. 8. Mitchell M. An Introduction to Genetic Algorithms. Cambridge, MA: The MIT Press; 1996. 9. Keedwell E, Narayanan A. Discovering gene networks with a neural-genetic hybrid. IEEE/ACM Trans Comput Biol Bioinform. 2005;2:231–242. 10. Lee WP, Tzou WS. Computational methods for discovering gene networks from expression data. Brief Bioinform. 2009;10:408–423. 11. Gonzalez-Álvarez DL, Vega-Rodríguez MA, Gómez-Pulido JA, Sanchez- Pérez JM. Predicting DNA motifs by using evolutionary multiobjective optimization. IEEE T Syst Man Cy C. 2012;42:913–925. 12. Fan Y, Wu W, Yang J, Liu R. An algorithm for motif discovery with iteration on lengths of motifs. IEEE/ACM Trans Comput Biol Bioinform. 2015;12:136–141. 13. Chaves-González JM, Vega-Rodríguez MA. DNA strand generation for DNA computing by using a multi-objective differential evolution algorithm. Biosystems. 2014;116:49–64. 14. Alba E, Troya JM. A survey of parallel distributed genetic algorithms. Complexity. 1999;4:31–52. 15. Cantú-Paz E. Efficient and Accurate Parallel Genetic Algorithms. Dordrecht, The Netherlands: Kluwer Academic Publishers; 2000. 16. Gong YJ, Chen WN, Zhan ZH, et al. Distributed evolutionary algorithms and their models: a survey of the state-of-the-art. Appl Soft Comput. 2015;34:286–300. 17. González DL, De Vega FF. On the intrinsic fault-tolerance nature of Parallel Genetic Programming. In: Proceedings of the 15th EUROMICRO International Conference on Parallel, Distributed and Network-based Processing; February 7-9, 2007:450–456; Naples. New York, NY: IEEE. 18. Hidalgo JI, Lanchares J, Fernández de Vega F, et al. Is the island model fault tolerant? In: Proceedings of the 2007 GECCO Conference Companion on Genetic and Evolutionary Computation; July 7-11, 2007:27–37; London, England. New York, NY: ACM. 19. Luo Z, Liu H. Cellular genetic algorithms and local search for 3-SAT problem on graphic hardware. In: Proceedings of the IEEE International Conference on Evolutionary Computation; July 16-21, 2006:2988–2992; Vancouver, BC, Canada. New York, NY: IEEE. 20. Li JM, Wang XJ, He RS, Chi Z-X. An efficient fine-grained parallel genetic algorithm based on GPU-accelerated. In: Proceedings of the IFIP International Conference on Network and Parallel Computing Workshops; September 18-21, 2007:857–864; Liaoning, China. New York, NY: IEEE. 21. Luong TV, Melab N, Talbi EG. GPU-based island model for evolutionary algorithms. In: Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation; July 7-11, 2010:1089–1096; Portland, OR. New York, NY: IEEE. 22. Ben-Shalom R, Aviv A, Razon B, Korngreen A. Optimizing ion channel models using a parallel genetic algorithm on graphical processors. J Neurosci Methods. 2012;206:183–194. 23. Arenas MG, Romero G, Mora AM, Castillo PA, Merelo JJ. GPU Parallel Computation in Bioinspired Algorithms: A Review. Berlin, Germany: Springer; 2012: 113–134. 24. Jin Y, Hallinan J. Evolving gene regulatory networks. Biosystems. 2009;98: 2008–2009. 25. Esmaeili A, Jacob C. A multi-objective differential evolutionary approach toward more stable gene regulatory networks. Biosystems. 2009;98:127–136. 26. Nicolau M, Schoenauer M. On the evolution of scale-free topologies with a gene regulatory network model. Biosystems. 2009;98:137–148. 27. Noman N, Monjo T, Moscato P, Iba H. Evolving robust gene regulatory networks. PLoS ONE. 2015;10:1–22. 28. Eggenberger P. Evolving morphologies of simulated 3d organisms based on differential gene expression. In: Fourth European Conference on Artificial Life, 28–31 July 1997:205–213. Brighton, UK: The MIT Press. 29. Banzhaf W. On the dynamics of an artificial regulatory network. In: Advances in Artificial Life: European Conference on Artificial Life; September 14-17, 2003:217–227; Dortmund, Germany. New York, NY: Springer. 30. Knabe J, Nehaniv CL, Schilstra M, et al. Evolving biological clocks using genetic regulatory networks. In: Artificial Life X: Proceedings of the Tenth International Conference on the Simulation and Synthesis of Living Systems; 3–6 June 2006. Bloomington, IN: The MIT Press. 31. Guo H, Meng Y, Jin Y. A cellular mechanism for multi-robot construction via evolutionary multi-objective optimization of a gene regulatory network. Biosystems. 2009;98:193–203. 32. Joachimczak M, Wróbel B. Processing signals with evolving artificial gene regulatory networks. In: Artificial Life XII: Proceedings of the Twelfth International Conference on the Synthesis and Simulation of Living Systems; 19–23 August 2010:203–212. Odense, Denmark. 33. Cussat-Blanc S, Harrington K, Pollack J. Gene regulatory network evolution through augmenting topologies. IEEE T Evolut Comput. 2015;19:823–837. 34. de Jong H. Modeling and simulation of genetic regulatory systems: a literature review. J Comput Biol. 2002;9:67–103. 35. Leclerc RD. Survival of the sparsest: robust gene networks are parsimonious. Mol Syst Biol. 2008;4:213. 36. Keckler SW, Dally WJ, Khailany B, Garland M, Glasco D. GPUs and the future of parallel computing. IEEE Micro. 2011;31:7–17. 37. Strohmaier E, Dongarra J, Simon H, et al. TOP500 supercomputing sites. http://www.top500.org/ 38. Halfhill TR. Parallel processing with CUDA. Microprocessor Report. 28 January 2008:1–8. 39. Nvidia. CUDA toolkit documentation. Technical Report. http://docs.nvidia.com/ cuda/. Published 2015. 40. Luque G, Alba E. Parallel Genetic Algorithms: Theory and Real World Applications (Vol 367). Berlin, Germany: Springer; 2011. 41. Krömer P, Platoš J, Snášel V. Nature-inspired meta-heuristics on modern GPUs: state of the art and brief survey of selected algorithms. Int J Parallel Prog. 2014;42:681–709. 42. Cecilia JM, Nisbet A, Amos M, García JM, Ujaldón M. Enhancing GPU parallelism in nature-inspired algorithms. J Supercomput. 2013;63:773–789. 43. Wahib M, Munawar A, Munetomo M, Akama K. Optimization of parallel Genetic Algorithms for nVidia GPUs. In: Proceedings of the 2011 IEEE Congress of Evolutionary Computation (CEC); June 5-8, 2011:803–811; New Orleans, LA. New York, NY: IEEE. 44. Sonnenburg S, Braun ML, Ong CS, et al. The need for open source software in machine learning. J Mach Learn Res. 2007;8:2443–2466. 45. Nvidia. NVIDIA visual profiler. https://developer.nvidia.com/nvidia-visual- profiler. Published 2016. 46. Jaros J. Multi-GPU island-based genetic algorithm for solving the knapsack problem. In: Proceedings of the 2012 IEEE Congress on Evolutionary Computation; June 10-15, 2012:1–8; Brisbane, QLD, Australia. New York, NY: IEEE. 47. Nvidia. NVIDIA thrust library. https://developer.nvidia.com/thrust. Published 2016. 48. Nvidia. CUDA toolkit 4.2 CURAND guide. Technical Report. https://developer. download.nvidia.com/compute/DevZone/docs/html/CUDALibraries/doc/ CURAND_Library.pdf. Published 2012. 49. Lässig J, Sudholt D. Design and analysis of migration in parallel evolutionary algorithms. Soft Comput. 2013;17:1121–1144. 50. Muelas S, Peña JM, Robles V, et al. Machine learning to analyze migration parameters in parallel genetic algorithms. Adv Soft Comp. 2007;44:199–206.