6.3.1 Mass-Action-Based Modeling of Gene Regulation

Given a biochemical system η consisting of some molecular species (e.g., genes, proteins, metabolites) and their mutual interactions, a mechanistic reaction-based model (RBM) of η can be formally defined by specifying the set of different molecular species occurring in η and the set of biochemical reactions. A reaction , μ = 1, …, M, is usually given in the form R_μ: a₁S₁ + ⋅⋅⋅ + a_NS_N → b₁S₁ + ⋅⋅⋅ + b_NS_N, where are called the stoichiometric coefficients of R_μ. The species occurring on the left-hand (right-hand) side of R_μ are called reagents (products, respectively). A numerical value is associated with R_μ, representing the rate constant of that reaction.

By considering the values of the stoichiometric coefficients of the reactions in , we can identify some particular reactions:

• if the stoichiometric coefficients a_i of the reagents are null for all i = 1, …, N, then R_μ is called a source reaction and it will be denoted as ∅ → products. Similarly, if the stoichiometric coefficients b_i of the products are null for all i = 1, …, N, then R_μ is called a sink (or degradation) reaction and it will be denoted as reagents → ∅. Reactions of the form ∅ → ∅, where a_i = 0 and b_i = 0 for all i = 1, …, N, are not considered in RBMs;
• meaningless reactions of the form a_iS_i → b_iS_i are excluded from any RBM since they correspond to unfeasible biochemical processes where a_i molecules of species S_i are converted into b_i molecules of the same species;
• at most second-order reactions (i.e., chemical reactions where at most two molecules occur as reactants) are considered in RBMs since third-order (or any higher order) reactions have a probability to occur almost equal to zero, as they would require the simultaneous collision of three (or more) reactant molecules.

In this chapter, we are interested, in particular, in reaction-based gene regulation models (GRMs), which formally describe the biochemical reactions involved in gene transcription and translation. To properly define these reactions, we assume that the set of molecular species is given by the union of two disjoint sets of species, which we denote as and , for some . We assume that |Γ| |Σ|, where | · | represents the cardinality of the set.

The set Γ represents the genes that are available for the synthetic engineering of gene regulation circuits. We assume that each element γ_i in Γ is strictly associated with a messenger RNA (mRNA), the product of gene transcription, whose expression can be experimentally evaluated through cutting-edge technologies (e.g., microarray [20], qRT-PCR [34]). Thus, the cardinality of Γ is equivalent to the number of target time series used in cuGENED. By abuse of notation, we will identify γ_i either with the “gene” or with its “transcription product”, according to the context.

The set Σ represents a set of generic species, as well as their mutual chemical complexes, that are related to the processes of gene expression. These species are the actual effectors of gene regulation, and need to be included in any GRM to evaluate their influence on the emerging behavior of the circuit. An element σ_i ∈ Σ might represent any kind of gene product (e.g., protein) related to γ_i, or the molecular complex formed by the interaction of the gene product (acting as promoting/inhibiting transcription factor) with another gene γ_j ∈ Γ.

To better clarify the meaning of the sets Γ and Σ, we illustrate some types of reactions that formally describe the processes of gene regulation and that might be present in a GRM, along with other generic reactions mentioned above.

Reactions that represent gene expression (transcription/translation) can be written in the form γ_i → γ_i + σ_i. Given a gene γ_j ∈ Γ, its regulation by means of gene γ_i can be described either as γ_i + γ_j → σ_k or as σ_i + γ_j → σ_k, with i not necessarily distinct from j. Here, the species σ_k might represent two different elements: (1) the product of the regulation of γ_j (e.g., the protein encoded by γ_j); in this case we can simply say that σ_k is equal to σ_j; (2) the chemical complex between gene γ_j and its own regulator; in this case σ_k is a compact formalization for the species γ_iγ_j or σ_iγ_j. The inverse reactions of the type σ_k → σ_i + γ_j can then be used to describe the dissociation of the regulator σ_i from gene γ_j. The effective expression of a regulated gene γ_j—that takes place after the occurrence of some reactions of type σ_i + γ_j → σ_k—can be synthetically written in the form σ_k → σ_j. This reaction states that species σ_j, related to gene γ_j, is derived from species σ_k, which represents an intermediate molecular complex having γ_j as element (e.g., σ_k = γ_iγ_j or σ_k = σ_iγ_j, as mentioned above). Finally, reactions of the form σ_h → σ_g can also be used to represent the degradation, or any other generic transformation reaction (e.g., post-translational modification of a protein) of the species occurring in Σ.

The RE process of a given biological system consists in the automatic identification of the network of interactions among the molecular components of that system. In this task, a set of experimental time series measurements of some species occurring in the system can be exploited as target to drive the evolutionary inference of the network [51]. On the contrary, in this chapter, we assume that only an expected dynamical behavior of the network is given as target. Quite obviously, we assume that no specific laboratory measurements of the temporal evolution of chemicals that will constitute the circuit can be given before the circuit itself has been engineered. For this reason, in what follows, we do not formally specify the type and the unit measurements of the amount of chemicals, but we use the general term “level”. The level can indicate either the number of molecules or the concentration values that might then be actually measured by ad hoc laboratory methodologies after the construction of the circuit. For instance, if qRT-PCR is exploited to validate the functioning of the synthetically engineered circuit, then the level of the species in Γ represents the concentrations of the corresponding mRNA, measured by relative quantification with respect to the amount of some control housekeeping gene [37].

According to the law of mass-action¹, given a set of biochemical reactions —each one characterized by its own kinetic constant k_μ² —for each molecular species appearing either as reagent or product in some reaction in , it is possible to derive a rate equation that describes the variation in its concentration with respect to time. In other words, any given RBM or GRM can be formalized as an equivalent system of coupled (nonlinear) first-order ODEs [71]. This formalization can be exploited to determine the temporal evolution of the network (that is, its dynamics) by means of different numerical integration methods [7]. These algorithms require as input the set of ODEs, along with the set of kinetic constants and the initial concentrations of the chemical species. To this aim, in this chapter, we exploit cupSODA [48, 50], a GPU-powered tool that automatically derives a set of ODEs from RBMs defined according to the mass-action kinetics, and then exploits the numerical integration algorithm LSODA [59] to perform the deterministic simulation of the system dynamics. In Section 6.4, we show how to exploit CGP to carry out the ED of a GRM, which is then automatically converted to an ODEs system and simulated by means of cupSODA, considering the kinetic constants that are inferred by PSO.

6.3.2 Cartesian Genetic Programming

CGP is an EC method whose individuals are described by means of indexed graphs, rather than trees, having (sequentially numbered) nodes arranged in a Cartesian coordinate grid [43]. The genotype of each individual is represented by a sequence of node connections and functions in the grid, called Cartesian program (CP).

Formally, a CP is a 9-tuple {G, n_i, n_o, F, FN, n_r, n_c, n_n, l} where

• G is the genotype, encoded as fixed-length vector of integer numbers representing the connections from the input nodes to the outputs nodes of the grid;
• are the number of input and output nodes, respectively;
• F is a finite set of functions (for instance, elementary arithmetic operations as { +, −, ×, /});
• FN is a grid of functional nodes, sequentially indexed by rows and columns, each one containing a function from F;
• are the number of functional nodes appearing in each row and in each column of the grid, respectively, so that |FN| = n_rn_c;
• is the number of input connections in each functional node;
• is the so-called “levels back” parameter, a measure of the CP inter-connectivity which determines how many preceding columns (in the grid of functional nodes) can have their output connected to the functional nodes appearing in any given column of the grid.

The connections of a CP start from the input nodes and pass through the functional nodes, each one having a fixed number n_n of input connections. The length of the genotype G is equal to n_rn_c(n_n + 1) + n_o, that is, an integer number is assigned to each input connection and to the output connection of every functional node, as well as to each output node of the grid. The inter-connectivity of nodes generally exploited in CGP, and specifically considered in this chapter, is strictly feed-forward, meaning that nodes belonging to the same column of the grid cannot be connected to each other. In addition, any node can be either connected or disconnected; disconnected nodes represent non-coding genes in the genotype and are ignored in the phenotype. The phenotype of a given CP is given by the actual graph that the CP represents. In this chapter, the semantic of the phenotype corresponds to the set of reactions derived from the graph, that will constitute the GRM. An example of CP genotype and phenotype, in the context of GRMs, is given in Section 6.4.

In a CGP population, after the fitness function of each CP has been evaluated, the best candidate solutions are selected to generate the offspring by means of a mutation operator. During the iterative process, CGP evolves a population of individuals characterized by a set of expressions, each one formed by the composition of input nodes and functional nodes. The total number of possible expressions represented by a given CP is upper bounded by the number n_o of output nodes. If an output node is linked to a disconnected functional node, some internal nodes will not be part of any path connecting an input node to that output node, a circumstance that leads to the possible existence of different genotypes mapping to the same phenotype. These non-coding regions are important in CGP for three reasons: (1) they reduce the size of the phenotype, as stated above; (2) they reduce the effect of mutations since a mutation acting on a disconnected functional node will not contribute to the variation of the phenotype; (3) even if they are unused in the current CP, a mutation could suddenly connect some disconnected node to the rest of the graph, resulting in a possibly relevant change in the phenotype.

In this chapter, the candidate solutions evolve by using a (1+λ) ES [5], as described in Ref. [43]: all individuals are evaluated and the best one is selected as a parent for the next generation. Then, λ offspring are produced by means of random mutations, that is, random modifications of the integers which constitute the genotype of the parent individual. The proportion of genes that are mutated is determined by the mutation rate parameter ρ ∈ (0, 1). The ES methodology does not exploit any crossover mechanism.

In Section 6.4, we show how to exploit CGP to perform the ED of GRMs. In order to evaluate the GRM that each CP represents, we need to define a proper fitness function. As target of the optimization process, we exploit a set of artificial temporal data that represent the desired behavior of the synthetic gene regulation circuit. These artificial data informally correspond to the concentrations of the molecular species that might be experimentally measured after the construction of the circuit itself. Therefore, the fitness function is based on the comparison of the target series against a simulated dynamics of the model, which is generated by exploiting the cupSODA tool [48, 50] explained in Section 6.3.4. Since the GRM derived with CGP is not complete until a proper kinetic parameterization is given, the fitness evaluation of each candidate solution embeds a PE phase, which is performed by means of PSO.

6.3.3 Particle Swarm Optimization

PSO is a population-based meta-heuristic inspired by the collective movement of flock of birds and school of fish [29]. In PSO, a swarm of P candidate solutions (named particles) moves inside an M-dimensional Euclidean search space, cooperating in the identification of the optimal solution, whose quality is measured using a specified fitness function. The movement of the swarm is bounded inside a subspace of : if a particle is driven outside the boundaries of this subspace as a consequence of its movement, then it is relocated inside the feasible region by means of the damping mechanism described in Ref. [72].

The i-th particle of the swarm is defined by two vectors: the position of the particle and its velocity . During the iterative optimization phase, the velocity is updated as a result of two attractors: the best position ever found by the particle itself, and the best position identified by the whole swarm. These attractors are balanced by means of two settings specified by the user: the cognitive factor and the social factor . The components of both attractors are further modulated by two vectors R₁ and R₂ of random numbers uniformly sampled in [0, 1), which allows the particles to move stochastically and to avoid the entrapment into local minima. Moreover, to reduce the possibility of chaotic movements of particles, the change of velocity is weighted by an inertia factor . The final definition of the velocity update formula for the i-th particle is:

(6.1)

where ○ represents the component-wise multiplication operator between vectors. Once the velocity is updated, the new position of each particle is determined by calculating φ_i = φ_i + v_i.

In this chapter, a particle corresponds to a candidate kinetic parameterization of each GRM determined by CGP. Therefore, M is equal to the number of connected outputs in the CP, that is, the total number of reactions in the GRM.

6.3.4 General-Purpose GPU Computing

The Compute Unified Device Architecture (CUDA) is a parallel computing platform and programming model introduced by Nvidia in 2006, consisting in a framework suitable to exploit GPUs in general-purpose computational tasks (GPGPU computing). GPGPU computing represents a valuable alternative to the traditional high-performance computing infrastructures (e.g., clusters of machines) since GPUs are characterized by low costs and a reduced energy consumption, allowing the access to tera-scale computing on common workstations of mid-range price. Nevertheless, a direct porting of sequential code on the GPU is most of the times unfeasible, due to the innovative architecture and the intrinsic limitations of this technology. Therefore, a fully exploitation of the computational power and massive parallelism of the GPU is usually challenging [17].

CUDA combines the single instruction multiple data (SIMD) architecture with multi-threading, and automatically handles the conditional divergence between threads. However, this flexibility has a drawback since any divergence of the execution flow among threads results in a serialization of the execution, affecting the overall performances. Following the naming conventions used in CUDA, a C/C++ function, called kernel, is loaded from the host (the CPU) to the devices (one or more GPUs) and replicated in many copies named threads. Threads can be organized in three-dimensional structures named blocks which, in turn, are contained in three-dimensional grids (as schematized in Figure 6.1, left side). Whenever the host computer runs a kernel function, the GPU creates the corresponding grid and automatically schedules each block of threads on an available streaming multi-processor of the GPU, thus allowing a transparent scaling of performances on different devices.

GPUs are equipped with different types of memory. As described in Figure 6.1, right side, the GPU memory hierarchy consists in the global memory (accessible from all threads), the shared memory (accessible from threads of the same block), the local memory (registers and arrays, accessible from owner thread), and the constant memory (cached and not modifiable). The best performances in the execution of CUDA code can be achieved by exploiting the shared memory as much as possible. Unfortunately, the shared memory is a very limited resource (49,152 bytes for each multi-processor, since the introduction of the Fermi architecture) that brings about restrictions on the blocks’ size. On the contrary, the global memory is very large (thousands of megabytes) but suffers from high latencies. In order to mitigate this issue, starting from the Fermi architecture, the global memory has been equipped with a L2 cache (see Figure 6.2). Moreover, with the Fermi architecture, the programmer can balance 64 KB of fast on-chip memory between the shared memory and L1 cache, specifying two different configurations: 48 KB for the shared memory and 16 KB for L1 cache, or 16 KB for the shared memory and 48 for L1 cache. In addition, using the Kepler architecture (exploited in our work), a third and perfectly balanced configuration can be specified by assigning the same amount of memory (32 KB) both to the shared memory and L1 cache.

Despite the remarkable advantages concerning the computational speedup, computing with GPUs usually requires the development and the implementation of ad hoc algorithms since GPU-based programming substantially differs from CPU-based computing. As a consequence, scientific applications of GPUs might undergo the risk of remaining a niche for few specialists. To avoid such limitations, several packages and software tools were recently released (see, e.g., Refs. [24, 50, 69, 73]), so that users with no knowledge of GPUs hardware and programming can also access the high-performance computing power of graphics engines.

In this chapter, the simulations needed to execute the ED of GRMs were carried out by means of cupSODA, a GPU-powered simulator for biological systems that allows the efficient execution of high numbers of parallel deterministic simulations [50, 48]. cupSODA relies on a C version of the numerical integration algorithm LSODA [59], ported and adapted to the CUDA architecture. This tool requires an Nvidia video card and can be executed on Microsoft Windows, Linux, and Apple OS X-based operating systems. cupSODA was designed to the purpose of being a black-box simulator, that can be easily used without any programming skills. In particular, cupSODA automatically converts any mass-action-based RBM of a biological system into the corresponding system of ODEs.

cupSODA exploits the massive parallelism of CUDA architecture for the parallel execution of independent simulations—each one characterized by a different initial parameterization of the model—in each thread. cupSODA was designed to speed up the time-consuming tasks typical of computational biology [2, 10], which rely on the repetition of large numbers of simulations in perturbed conditions, generally realized by varying the initial concentrations of chemical species or the value of the kinetic constants. cupSODA takes advantage of the memory hierarchy of the GPU, by allocating the most frequently updated data (e.g., the state of the system) in the shared memory, and the unvarying data (e.g., the set of reactions in the RBM) in the constant memory.

cupSODA integrates a specific functionality for the evaluation of the fitness function since it allows to compare the outcome of the simulations with any available target data (e.g., the levels of chemical species at chosen time instants t₀, …, t_F). This is useful, for instance, to perform a PE task. To this aim, cupSODA invokes the LSODA kernel F − 1 times: each time the kernel is run over a (simulated) time interval of length Δt = t_i − t_{i − 1}, i = 1, …, F, and the concentration values of the output species are stored at the end of each Δt. Once the concentrations are stored, cupSODA provides a set of metrics (e.g., Equation 6.2) that can be exploited to evaluate in parallel the fitness function of each particle.

CGP phase.

In cuGENED, CGP is implemented using F = { +, −} as the set of functions. The composition of input nodes and functional nodes allows CGP to encode complex expressions for each genotype G. The expressions can be then converted, by means of symbolic manipulation, into arithmetical equations as described hereafter. In the context of ED of GRMs, the set of input nodes is determined by the species occurring in Γ and Σ, so that n_i = |Γ| + |Σ|. Since the exact regulatory interactions necessary to reproduce the desired behavior of the synthetic circuit are to be determined, it follows that the number of chemical species belonging to the set Σ is unknown as well. Therefore, also the cardinality of Σ must be inferred by means of some heuristics. As a matter of fact, the value |Σ| has a relevant impact on the ED process since it determines the space of the possible GRMs that CGP can explore. If |Σ| is too small, then a circuit that is able to achieve a perfect fit of the desired dynamics might be impossible to be designed. On the contrary, if |Σ| is too large, then the evolutionary algorithm can take a longer time to converge to an optimal solution, or might suffer from overfitting. Nevertheless, in the second case, CGP is able to automatically exclude the unnecessary chemical species from the set of candidate solutions if they have no impact on the fitness value. In this work, we always use the value |Σ| = 4 · |Γ| as heuristic, because we assume that, for each target species in Γ, four intermediates are enough to model gene regulation mechanisms.

The number of input connections for each functional node is n_n = 2, because we only exploit addition and subtraction as functions for the construction of the expressions. The set of output nodes corresponds to the biochemical reactions that will be part of the inferred GRM. In these reactions, the species corresponding to the input nodes appear either as reactants, if a functional node containing function + is crossed, or as products, if a functional node containing function is crossed. In so doing, each expression encodes a single reaction of the GRM.

Generally, in CGP, the number of output nodes n_o is fixed. In the context of GRMs, this implicitly means that the actual number of reactions should be known before the optimization takes place, which is clearly unreasonable. Since the GRM itself is the goal of the optimization, in this chapter, the value n_o represents an upper bound to the number of reactions that will appear in the model. Hence, a proper choice for n_o is fundamental: it should be large enough to include all the necessary regulatory mechanisms, and small enough to avoid bloating and overfitting. As heuristic for the selection of n_o, we suggest to fix its value equal to the number of all possible reactions that can involve the set of genes and the set of generic chemical species, that is, n_o = |Γ| · |Σ| + |Γ| + |Σ|. Furthermore, in this implementation, we set l = n_c, meaning that output nodes can be connected to any functional or input node, and each functional node belonging to column j can be connected only to nodes between column j − 1 and the input nodes.

An example of CP related to the ED problem of GRMs is explained in Example 6.1 and graphically represented in Figure 6.3.

■ EXAMPLE 6.1

There exists a strict correspondence between the genotype of the CP and the connectivity of each node in the Cartesian coordinate grid. Functions + and − in the grid are represented by 0 and 1 in the genotype, respectively (the first number in each of the four triplets). Since + and − are binary operators, we use n_n = 2 input connections to each functional node. The grid of functional nodes consists in |FN| = n_rn_c = 4 nodes since the grid is composed of two rows and two columns. The genotype is generated by randomly drawing values related to the chemical species (0 for γ₀, 1 for σ₀, 2 for σ₁) and to the output of each functional node (3, … , 6). Then, these values are assigned to the input of the functional nodes and to the output nodes (numerical values with asterisk). Note that, in the construction of the genotype, given a functional node, only the values of previous input and functional nodes can be randomly selected and assigned to its input connections since l = n_c. The grid is parsed backwards, from the n_o = 2 output nodes R₁, R₂ to the n_i = 3 input nodes γ₀, σ₀, σ₁. The resulting equations are then automatically converted into an equivalent set of biochemical reactions. The grey functional nodes are involved in the equations since a path from the input nodes to the output nodes exists here. The white functional nodes, instead, are not connected in the grid and do not participate in any equation in this CP. Thus, their corresponding genes (gray numbers) are non-coding sequences.

The pseudo-code of our ED methodology based on CGP and PSO is reported in Algorithm 6.1. The evolutionary process begins with the creation of a population of I = 1 + λ random CPs (lines 2–7). The population evolves by means of an ES process which can be decomposed in three main steps:

1. in Step 1, the representation of each CP is converted into a GRM. Then, PSO is used to estimate the values of the kinetic constants and to assess the fitness value of the parameterized candidate network (lines 10–14);
2. in Step 2, the CPs are ranked according to the fitness values, to identify the best CP in the population (line 15);
3. in Step 3, a brand new population is formed by considering the best CP, together with the I − 1 offspring created by applying the mutation operator on the best CP (lines 16–30).

Figure 6.3 schematizes the conversion of a CP into the corresponding GRM (Step 1). The connections of each output node are followed backwards by recursively passing through the functional nodes. This process yields a derivation tree that is translated into an arithmetical equation, where positive terms (respectively, negative) are considered as reactants (products). Each equation produces a single candidate chemical reaction for the network. A candidate GRM η for the synthetic circuit is then obtained by repeating the algorithm for all output nodes.

The fitness values of all candidate solutions are calculated at the end of the PE process by means of PSO, so that a ranking of the CPs can be assessed and the best solution (best_CP) in the CGP population is identified (Step 2).

During each generation of the CGP, the new offspring solutions are obtained by applying a mutation operator to best_CP. The GRM η is compared to the rest of the population: if it is identical to some other CP, it is mutated again in order to achieve a heterogeneous population, and the process is repeated. In this chapter, we do not verify that η consists of a single connected component since not all chemical species in Σ necessarily need to be present in the GRM. The rationale behind this choice is that, in the ED methodology presented here, we do not specify a priori the exact number of chemical species that should occur in the system (we only provide an upper bound). Therefore, in this methodology also the number of chemical species in Σ undergoes the optimization process. This approach is different from typical methods used for the RE problem. In particular, it differs from the RE method based on CGP and PSO that we previously proposed in Ref. [51]: in that case, the number of chemical species participating in the system was known a priori; therefore, it was mandatory for the candidate solutions to consist in a single connected component.

When a new population of CGP is generated, for each candidate solution we perform a consistency check (lines 24–28 in Algorithm 6.1). Specifically, we verify that each reaction in the GRM η obeys the conditions described in Section 6.3.1, together with the additional condition that no identical reactions appear in η. In particular, since we consider at most second-order reactions, we can limit the length of CP expressions and reduce the bloating by setting n_c 3. In general, if a reaction R is not consistent, we remove it from the network, so that η = η∖{R}. Eventually, when a new consistent CP is produced, the individual is inserted into the new population.

To improve the convergence speed, at the end of Step 3 we apply the elitist selection by adding a non-mutated copy of the best candidate network found so far. This solution becomes the (I + 1)th individual of the population (lines 31–33 in Algorithm 6.1). When the new population is formed, we calculate the fitness value of each individual by executing a PE by means of PSO, as described hereby. This iterative process is repeated for generations. Finally, the GRM with the best fitness, along with its parameterization, is chosen as the result of the ED problem.

PSO phase.

For the execution of the PE task, each particle in the PSO corresponds to a candidate kinetic parameterization of a GRM inferred by CGP. Namely, for each candidate GRM, we exploit a whole swarm of particles to determine the reaction constants that best fit the desired target behavior for that GRM. To this aim, the PSO adopts the following settings: P = 64 particles, randomly generated using a logarithmic distribution to better distribute the values of kinetic constants over different orders of magnitude; c_soc = c_cog = 2.05, as suggested in Ref. [13]; inertia linearly decremented from w = 0.9 to w = 0.4; velocity vector automatically clamped to a maximum intensity, equal to 10% of the maximum point-to-point distance in the search space; search space bounded between 10^{− 5} and 10² (for each kinetic parameter) with damping boundary conditions [72]. The choice of the upper bound of the search space for the kinetic constants is fundamental for the entire optimization process since a limited range of variation might exclude the global optimum. However, to perform PE of biochemical systems consisting in a small set of reactions, and characterized by molecular species with low concentrations, limiting the search space of parameters is a good practice for the following reason. Very high values of kinetic constants could lead to undesired dynamics in which chemical species are entirely consumed in the first time instants of the simulation.

The movement of the swarm during the PE phase is determined according to the fitness function. In this chapter, for each particle φ_i, i = 1, …, P, which encodes the kinetic parameterization of a GRM η, the fitness is defined as the normalized distance between the values of the simulated dynamics of gene expression levels in η and the desired target dynamics. Formally, given the set Γ of genes that are available for the synthetic engineering of the gene regulation circuit and the set of reactions of the GRM inferred by CGP, we denote by X_h(t_c) the expected measurement at time t_c of the hth chemical species in Γ and by its simulated value sampled at time t_c. The value is obtained by means of cupSODA, using the kinetic constants contained in particle φ_i. The fitness function of particle φ_i is therefore given by:

(6.2)

where corresponds to the number of time instants that are arbitrarily sampled in the target dynamics of the GRM.

In this chapter, the PSO algorithm is halted after iterations. Then, for each GRM η, the value —that is, the fitness of the best solution g = (k₁, …, k_M) determined within the swarm of all candidate parameterizations of η—is taken as the fitness of the inferred GRM η. The GRM characterized by the minimum fitness value among all GRMs inferred by CGP is eventually chosen as the best solution of the whole ED problem.

GPU implementation.

Fitness evaluations are computationally expensive since multiple simulations (one for each particle) must be performed during all iterations of the PE on each candidate network inferred by CGP. More precisely, it takes fitness evaluations to perform a whole ED process. Anyway, all simulations during each PE phase are mutually independent and can be straightforwardly accelerated by means of a parallel architecture. In order to take advantage of the parallelism of modern GPUs, we exploit the GPU-accelerated simulator cupSODA [48, 50] and launch P threads, which perform the simulations and calculate the fitness functions defined in Equation 6.2. This way, the impact of fitness evaluations on the overall running time decreases to .

The rest of the methodology (CGP and PSO) is implemented using the Python language (version 2.7) and is executed in a strictly sequential fashion. PSO invokes cupSODA to assess the fitness values by means of synchronous subprocess calls.

6.5 RESULTS

To test the feasibility and the effectiveness of cuGENED, we performed the ED of synthetic circuits composed by two and three genes, using as target of the optimization process a desired temporal dynamics of a small set of (in silico generated) mRNA levels. The target dynamics were created from scratch to be complex enough to require the ED of GRMs possibly characterized by multiple regulation mechanisms and intermediate species.

As previously mentioned, the results of the ED process can be only evaluated by considering the desired dynamics. No target GRM can be exploited to assess the quality of the inferred network since the goal of the ED process is to determine a synthetic circuit that still has to be constructed in laboratory. Therefore, to prove the effectiveness of cuGENED, only comparisons between the expected and the simulated dynamics of the inferred GRM are presented hereby.

All the following tests were performed on a workstation with a CPU Intel Core i7-2600, with a clock frequency of 3.4 GHz, and with a GPU Nvidia GeForce GTX 590, running OS Windows 7 64 bit. The settings of the LSODA integrator used for the simulation of the dynamics of the candidate solutions are: relative error equal to 1 · 10^{− 10}, absolute error equal to 1 · 10^{− 10}, maximum number of integration steps equal to 10,000.

6.5.1 ED of Synthetic Circuits with Two Genes

For the ED of synthetic circuits consisting of two genes, we set |Γ| = 2 and |Σ| = 8, so that the inferred GRM can contain at most n_i = 10 chemical species and n_o = 26 reactions (according to the heuristics described in Section 6.4). The desired target dynamics of the species in Γ are described by two sigmoidal curves, characterized by different slopes. This behavior corresponds to a down-regulation of both genes, whose constitutive expression decreases in time by the action of some mutual regulatory mechanism, that has to be inferred by cuGENED.

Some preliminary tests were executed to analyze the influence of cuGENED settings on the ED process. In the first test, we investigated the impact of the parameter n_r—that is, the number of rows in the grid of functional nodes—on the convergence speed (Figure 6.4). According to our results for the two-genes system under investigation, too small (n_r = 8) or too large (n_r = 14) values of n_r yield a worse convergence of cuGENED, while the best results are achieved with n_r = 10. Moreover, by analyzing the early phases of the CGP process we can observe that the higher the n_r value the faster the convergence. Nevertheless, it is the intermediate value n_r = 10 that achieves the best convergence speed, confirming a similar result presented in Ref. [51]. In addition, Figure 6.4 highlights a feature that is typical of methodologies embedding PSO into CGP, as previously discussed in Ref. [51]: due to the stochasticity of PSO, two different PE executions on the same network usually lead to two different kinetic parameterizations which, in turn, are likely to have two different fitness values, explaining the large fluctuations in the convergence speed plots.

As a second test, we analyzed the impact of the parameter ρ, which determines how many elements of a CP genotype are mutated during the offspring generation. Figure 6.5 shows the convergence of cuGENED with different values of ρ, highlighting that the best choice is ρ = 0.2. This setting represents the best trade-off between exploration and exploitation of the search space. Indeed, if the value of ρ is too small, CGP is not able to properly explore the search space and cannot converge to individuals characterized by a good fitness value. On the contrary, if the value of ρ is too large, the mutation easily disrupts the structure of the best individual, yielding random individuals with a large fitness value and reducing the advantage of an evolution-guided exploration of the search space.

The best setting identified for cuGENED (i.e., n_r = 10 and ρ = 0.2) was then exploited to derive the GRM η_2γ that achieves the best fit with the expected behavior of the two-genes system. Figure 6.6 shows the comparison between the desired target dynamics (dots) and the simulated dynamics (lines) of η_2γ. This network is represented in Figure 6.7 and consists of the following reactions:

where the numerical values above the arrows correspond to the best kinetic parameterization determined by PSO.

All the reactions in η_2γ are biologically consistent and contribute well to the definition of a network of plausible genetic interactions, which could be implemented by means of synthetic biology techniques. For instance, in η_2γ we can observe the presence of a multi-product reaction (R₆) [14] and alternative splicing reactions (R₄ and R₅) [40]. Furthermore, η_2γ is characterized by a low fitness value, which indicates a consistent similarity with the expected target dynamics obtained by using the kinetic constants (k₁, …, k₈) inferred by PSO (Figure 6.6). Finally, it is worth noting that this optimal solution does not exploit all the chemical species in the set Σ: the CGP avoided the bloating and evolved a well-fitting GRM that only contains 4 out of the 8 possible species in Σ.

6.5.2 ED of Synthetic Circuits with Three Genes

As a further test, we performed the ED on a system consisting of three genes, using the best settings of cuGENED identified during the preliminary tests (Section 6.5.1). For the ED of this synthetic circuit, we set |Γ| = 3 and |Σ| = 12, so that the inferred GRM can contain at most n_i = 15 chemical species and n_o = 51 reactions (according to the heuristics described in Section 6.4).

The desired target dynamics of the species in Γ are described by two sigmoidal curves, characterized by different slopes, and a monotonically increasing curve (Figure 6.8, dots). This behavior corresponds to a down-regulation of two genes and an up-regulation of a third gene, whose underlying regulatory interplay is to be inferred by means of cuGENED.

The best GRM η_3γ inferred by cuGENED is represented in Figure 6.9 and consists in the following reactions:

where the numerical values above the arrows correspond to the best kinetic parameterization determined by PSO.

The result in Figure 6.8 shows that η_3γ almost perfectly fits the desired dynamics. In particular, the target and simulated dynamics for γ₀ and γ₁ are perfectly overlapped, while the curve of γ₂ slightly diverges from the expected behavior. Nevertheless, the simulated behavior of γ₂ is qualitatively and quantitatively similar to the target behavior, which means that, from a biological standpoint, the ED goal is fully achieved.

The interaction network in Figure 6.9 has two important characteristics, highlighted by the dashed lines, which represent side effects of cuGENED:

• the chemical species σ₂, which is a product of reaction R₅, is not involved in the rest of the network. Hence, this species can be removed from η_3γ without any impact on the system dynamics. We highlight that the role of the chemical species σ₃ is different from σ₂, because it is fundamental for the degradation of γ₀ as a consequence of reaction R₁;
• the reactants of reaction R₇ are species σ₆ and σ₁₀: while the latter is produced by reaction R₃, the former is not produced by any other reaction, so that the level of σ₆ remains fixed to 0 during the simulations. The consequence is that reaction R₇ will never take place and hence it can be removed from η_3γ, along with σ₆. On the contrary, species σ₁₁, the product of reaction R₇, cannot be removed from the network, being involved also in reaction R₆.

6.6 DISCUSSION

cuGENED consists in the integration of two evolutionary algorithms which work at two different levels: CGP to infer the structure of a GRN, and PSO to estimate the kinetic parameters of the reactions involved in the network. To the best of our knowledge, cuGENED represents the first attempt to tackle the ED of GRMs by means of CGP combined with PSO. This methodology is also the first attempt to automatically model a GRN assuming 0-knowledge on the set of molecular species. In this context, there are some issues that are worth considering.

In cuGENED, the ED problem is formalized without any a priori information on the structure of the system, defined in terms of molecular interactions, except for a desired target dynamics of a small subset of species. On the one hand, this information is sufficient to properly reconstruct a network whose dynamics fits the target data (see, for instance, Figure 6.6). On the other hand, this information alone is usually not enough to discriminate between possibly different network topologies [11], which could all be able to reproduce the same behavior used as target of the optimization. This well known issue in ED and RE is named indistinguishability of equivalent networks [68]. To lessen this computational difficulty, a good practice should be the multiple execution of the ED process, followed by an analysis of the obtained solutions to assess their laboratory feasibility. This analysis could as well be performed by means of automatic algorithms, which exploit a knowledge base of feasible reactions specifically defined by synthetic biologists.

With respect to other existing methodologies for the inference of GRNs, an advantage of cuGENED is that, thanks to the use of CGP, it allows to exploit a simple representation of the candidate GRMs. The genotype of CGP, indeed, is a fixed-length vector of integer values, whose corresponding phenotype is a set of human-readable chemical reactions. Moreover, cuGENED requires the estimation of only M kinetic parameters associated with the reactions, where, in general, M ≪ (N² + 1)/2. On the contrary, other existing GP methodologies generally necessitate complex parameterizations. For instance, S-tree-based methods require the estimation of 2N(N + 1) kinetic parameters [9], thus making this method hardly exploitable for systems involving many genes. Moreover, thanks to the explicit limitation to second-order reactions in GRMs, cuGENED does not generate families of candidate solutions having an identical structure but a different stoichiometry, which still represents an open issue in GP-based methodologies [39].

Another relevant strength of cuGENED is the ED of synthetic circuits formally defined as GRMs, which, in contrast to most existing approaches [8], are not based on arbitrary kinetics functions. The GRNs evolved by CGP and PSO consist in mass-action reaction-based models, which describe the biological processes in terms of simple molecular interactions, and not as reaction rate approximations based on some kind of chemical assumption. As a matter of fact, Hill functions were shown not to represent physically realistic reaction schemes [70], especially for gene expression processes [63], while the equilibrium or quasi steady-state assumptions at the basis of Michaelis–Menten constants are valid only in specific conditions [46].

A current drawback of cuGENED is that the choice of its initial setting is relevant to converge to an optimal solution since both CGP and PSO are not settings-free algorithms. We relied on common literature settings for the PSO, while we analyzed the impact of different parameter values on the performances of CGP (i.e., the mutation rate ρ and the number of rows n_r). As a future improvement of cuGENED, we will investigate the relationship between these two values, the number of genes and the complexity of the desired target dynamics, in order to derive a heuristic for the automatic selection of the most appropriate settings.

cuGENED relies on the nested execution of two evolutionary algorithms. As a consequence, its overall computational cost can be relevant, being proportional to the number of generations of CGP, to the number of iterations of PSO, to their population size, to the number of reactions involved in the candidate networks, and to the time length of the simulations. The largest part of the running time is due to the simulation of the dynamics of the candidate solutions, which is a fundamental step for the fitness computation. However, the computational cost can be strongly reduced by exploiting the GPU-powered cupSODA simulator [50]. In so doing, we realized a parallel execution of the simulations and the fitness calculations of all candidate solutions, achieving a 24.5 × speedup with respect to a strictly sequential execution of the same tasks. It is worth noting that the computational time could be further reduced by executing parallel ED instances. To the best of our knowledge, cuGENED also represents the first attempt ever in exploiting GPUs to perform an accelerated ED of GRMs.

6.7 CONCLUSIONS AND FUTURE PERSPECTIVES

The idea of a circuit-like connectivity between biological parts was postulated for the first time in the 1960s [45]. This intuition lead to several attempts to properly formalize biological regulation systems through mathematical models [22, 23, 28, 64] and analyze the cellular pathways under investigation by exploiting electrical circuit analogies [41, 42]. This research field was favored by the progression of molecular biology and genomics, which provided several methods and the necessary knowledge to physically assemble biomolecular components. As a matter of fact, it is nowadays possible to properly evaluate the interaction among different genes through expression data coming from high-throughput techniques like microarrays [20] and real-time PCR [34]. It is also possible to engineer in vitro or in vivo customized signaling circuits [36]. The proof of concept that a computing-like behavior could be applied to biological systems was the design of the first synthetic gene networks, realized by using engineering-based methodologies [16, 19]. Lately, several network models were proposed, all of them integrating biochemical pathway information and expression data [1, 6, 38, 58].

In this context, we developed cuGENED, a computational methodology for the automatic design of GRNs characterized by a predefined dynamical behavior. In particular, cuGENED allows to derive a mechanistic model of a gene regulation system, modeled as a parameterized set of biochemical reactions. The reactions describe the processes related to gene expression, and they involve a set of molecular species (e.g., genes, mRNA) whose dynamics represent the target behavior of the optimization.

It is often the case that a gene can regulate either itself (autoregulation) or a single other gene. This is due to the fact that a generic regulatory species might be able to interact with a single component of the gene network through a unique DNA sequence. This means that all species could be extremely selective for their molecular target (i.e., each species can bind to a unique DNA sequence), so that no unknown interactions can occur between the regulators and their corresponding molecular targets [18]. Therefore, the molecular target of a regulator could be either the gene itself or another gene, for example, in a very selective feedback system. Thus, as future developments of this work, we plan to improve the ED methodology as follows:

• we will include the verification of the requirements about selective regulation to narrow the search space of candidate solutions, possibly obtaining GRMs that could be easier to implement with laboratory techniques;
• we will exploit the massive parallelism provided by GPUs to perform many parallel ED processes, in order to collect a set of optimal GRMs. This set can then be analyzed to identify the GRMs that show a better selective regulation (i.e., the GRMs that better fulfill the requirements stated above).

We also plan to introduce an automatic mechanism to detect the unessential species (e.g., products that have no impact on the system dynamics), as well as the reactions that cannot take place (i.e., reactions whose reactants have zero concentration, like those highlighted in Figure 6.9). In addition, we will consider the fact that the inference process might converge to a GRM whose simulation perfectly fits the desired behavior, though some reactions appearing in the GRM might not be biologically plausible. As mentioned in Section 6.6, a solution to this problem might come from the inclusion of some knowledge domain constraints in the ED, in order to automatically remove such reactions from the candidate GRMs.

Finally, in cuGENED we have not explicitly taken into account the biological noise that is typical of gene expression [15]. Actually, some species involved in GRNs can have very low intracellular amounts (in the order of tens or a few hundreds of molecules), so that a stochastic simulation methodology could be preferable with respect to deterministic simulations to correctly reproduce any noise-induced emergent phenomena. cuGENED can be straightforwardly extended to support also stochastic simulation algorithms for the simulation of the temporal dynamics of candidate solutions [21, 52, 62]. Indeed, as a further development, we will investigate the possibility of the automatic reconstruction of networks in presence of intrinsic noise, oscillations, or multistability phenomena.

REFERENCES

1. T. Akutsu, S. Miyano, and S. Kuhara. Identification of genetic networks from a small number of gene expression patterns under the Boolean network model. In Pacific Symposium on Biocomputing, volume 4, pages 17–28, 1999.
2. B. B. Aldridge, J. M. Burke, D. A. Lauffenburger, and P. K. Sorger. Physicochemical modelling of cell signalling pathways. Nature Cell Biology, 8:1195–1203, 2006.
3. S. Ando, E. Sakamoto, and H. Iba. Evolutionary modeling and inference of gene network. Information Sciences, 145:237–259, 2002.
4. D. Besozzi, P. Cazzaniga, G. Mauri, D. Pescini, and L. Vanneschi. A comparison of genetic algorithms and particle swarm optimization for parameter estimation in stochastic biochemical systems. In C. Pizzuti, M.D. Ritchie, and M. Giacobini, editors, Evolutionary Computation, Machine Learning and Data Mining in Bioinformatics (Proceedings EvoBIO 2009), volume 5483 of LNCS, pages 116–127. Springer, 2009.
5. H. Beyer and H. Schwefel. Evolution strategies—a comprehensive introduction. Natural Computing, 1:3–52, 2002.
6. J. M. Bower and H. Bolouri. Computational Modeling of Genetic and Biochemical Networks. MIT Press, 2001.
7. J. C. Butcher. Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, 2003.
8. H. Cao, L. Kang, Y. Chen, and J. Yu. Evolutionary modeling of systems of ordinary differential equations with genetic programming. Genetic Programming and Evolvable Machines, 1:309–337, 2000.
9. D. Cho, K. Cho, and B. Zhang. Identification of biochemical networks by S-tree based genetic programming. Bioinformatics, 22(13):1631–1640, 2006.
10. I. C. Chou and E. O. Voit. Recent developments in parameter estimation and structure identification of biochemical and genomic systems. Mathematical Biosciences, 219:57–83, 2009.
11. G. Craciun and C. Pantea. Identifiability of chemical reaction networks. Journal of Mathematical Chemistry, 44(1):244–259, 2008.
12. L. Demattè and D. Prandi. GPU computing for systems biology. Briefings in Bioinformatics, 11(3):323–333, 2010.
13. A. Dräger, M. Kronfeld, M. J. Ziller, J. Supper, H. Planatscher, and J. B. Magnus. Modeling metabolic networks in C. glutamicum: a comparison of rate laws in combination with various parameter optimization strategies. BMC Systems Biology, 3(5), 2009.
14. X. Du, J. Wang, H. Zhu, L. Rinaldo, K. Lamar, A. C. Palmenberg, C. Hansel, and C. M. Gomez. Second cistron in CACNA1A gene encodes a transcription factor mediating cerebellar development and SCA6. Cell, 154(1):118–133, 2013.
15. A. Eldar and M. B. Elowitz. Functional roles for noise in genetic circuits. Nature, 467(7312):167–173, 2010.
16. M. B. Elowitz and S. Leibler. A synthetic oscillatory network of transcriptional regulators. Nature, 403(6767):335–338, 2000.
17. R. M. Farber. Topical perspective on massive threading and parallelism. Journal of Molecular Graphics and Modelling, 30:82–89, 2011.
18. N. C. Garbett and J. B. Chaires. Binding: a polemic and rough guide. Methods in Cell Biology, 84:1–23, 2008.
19. T. S. Gardner, C. R. Cantor, and J. J. Collins. Construction of a genetic toggle switch in Escherichia coli. Nature, 403(6767):339–342, 2000.
20. G. Gibson. Microarray analysis. PLoS Biology, 1(1):e15, 2003.
21. D. T. Gillespie. Stochastic simulation of chemical kinetics. Annual Review of Physical Chemistry, 58:35–55, 2007.
22. L. Glass. Classification of biological networks by their qualitative dynamics. Journal of Theoretical Biology, 54(1):85–107, 1975.
23. L. Glass and S. A. Kauffman. The logical analysis of continuous, non-linear biochemical control networks. Journal of Theoretical Biology, 39(1):103–129, 1973.
24. M. J. Harvey and G. De Fabritiis. A survey of computational molecular science using graphics processing units. Wiley Interdisciplinary Reviews: Computational Molecular Science, 2(5):734–742, 2012.
25. S. Hoops, S. Sahle, R. Gauges, C. Lee, J. Pahle, N. Simus, M. Singhal, L. Xu, P. Mendes, and U. Kummer. COPASI—a COmplex PAthway SImulator. Bioinformatics, 22:3067–3074, 2006.
26. H. Iba. Inference of differential equation models by genetic programming. Information Sciences, 178:4453–4468, 2008.
27. E. Jones, T. Oliphant, and P. Peterson. SciPy: Open source scientific tools for Python, 2001.
28. S. A. Kauffman. The large scale structure and dynamics of gene control circuits: An ensemble approach. Journal of Theoretical Biology, 44(1):167–190, 1974.
29. J. Kennedy and R. C. Eberhart. Particle swarm optimization. In Proceedings of the IEEE International Conference on Neural Networks, volume 4, pages 1942–1948, Piscataway, NJ, 1995.
30. J. Kitagawa and H. Iba. Identifying metabolic pathways and gene regulation networks with evolutionary algorithms. In Gary B. Fogel and David W. Corne, editors, Evolutionary Computation in Bioinformatics, pages 255–278. San Francisco: Morgan Kaufman, 2003.
31. J. R. Koza. Genetic Programming: On the Programming of Computers by Means of Natural Selection. The MIT Press, 1992.
32. J. R. Koza, W. Mydlowec, G. Lanza, J. Yu, and M. A. Keane. Reverse engineering of metabolic pathways from observed data using genetic programming. In Pacific Symposium on Biocomputing, volume 6, pages 434–445, 2001.
33. J. R. Koza, W. Mydlowec, G. Lanza, J. Yu, and M. A. Keane. Reverse engineering of metabolic pathways from observed data using genetic programming. Technical Report SMI-2000-0851, Stanford University, Stanford, California (USA), 2001.
34. M. Kubista, J. M. Andrade, M. Bengtsson, A. Forootan, J. Jonák, K. Lind, R. Sindelka, R. Sjöback, B. Sjögreen, L. Strömbom, A. Ståhlberg, and N. Zoric. The real-time polymerase chain reaction. Molecular Aspects of Medicine, 27(2–3):95–125, 2006.
35. T. Lenser, T. Hinze, B. Ibrahim, and P. Dittrich. Towards evolutionary network reconstruction tools for Systems Biology. In E. Marchiori, J. H. Moore, and J. C. Rajapakse, editors, Evolutionary Computation, Machine Learning and Data Mining in Bioinformatics (Proceedings EvoBIO 2007), volume 4447 of LNCS, pages 132–142. Springer, 2007.
36. W.A. Lim. Designing customized cell signalling circuits. Nature Reviews Molecular Cell Biology, 11(6):393–403, 2010.
37. J. Logan, K. Edwards, and N. Saunders, editors. Real-Time PCR: Current Technology and Applications. Caister Academic Press, 2009.
38. Y. Maki, D. Tominaga, M. Okamoto, S. Watanabe, and Y. Eguchi. Development of a system for the inference of large scale genetic networks. In Pacific Symposium on Biocomputing, volume 6, pages 446–458, 2001.
39. D. Marco, C. Shankland, and D. Cairns. Evolving Bio-PEPA algebra models using Genetic Programming. In Proceedings of the 2012 Annual Conference on Genetic and Evolutionary Computation, pages 177–183. ACM, 2012.
40. A. J. Matlin, F. Clark, and C. W. J. Smith. Understanding alternative splicing: towards a cellular code. Nature Reviews Molecular Cell Biology, 6(5):386–398, 2005.
41. H. H. McAdams and A. Arkin. Gene regulation: Towards a circuit engineering discipline. Current Biology, 10(8):318–320, 2000.
42. H. H. McAdams and L. Shapiro. Circuit simulation of genetic networks. Science, 269(5224):650–656, 1995.
43. J. Miller and P. Thomson. Cartesian genetic programming. In R. Poli, W. Banzhaf, W. B. Langdon, J. Miller, P. Nordin, T. C. Fogarty, editors, Proceedings of the Third European Conference on Genetic Programming (EuroGP2000), volume 1802 of LNCS, pages 121–132. Springer, 2000.
44. C. G. Moles, P. Mendes, and J. R. Banga. Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome Research, 13:2467–2474, 2003.
45. J. Monod and F. Jacob. General conclusions: teleonomic mechanisms in cellular metabolism, growth, and differentiation. In Cold Spring Harbor Symposia on Quantitative Biology, volume 26, pages 389–401. Cold Spring Harbor Laboratory Press, 1961.
46. D. L. Nelson and M. M. Cox. Lehninger Principles of Biochemistry. W. H. Freeman Company, 2004.
47. M. S. Nobile, D. Besozzi, P. Cazzaniga, and G. Mauri. The foundation of evolutionary Petri nets. In G. Balbo and M. Heiner, editors, Proceedings of the 4th International Workshop on Biological Processes & Petri Nets (BioPPN 2013), volume 988, pages 60–74. CEUR Workshop Proceedings, 2013.
48. M. S. Nobile, D. Besozzi, P. Cazzaniga, and G. Mauri. GPU-accelerated simulations of mass-action kinetics models with cupSODA. The Journal of Supercomputing, 69(1):17–24, 2014.
49. M. S. Nobile, D. Besozzi, P. Cazzaniga, G. Mauri, and D. Pescini. A GPU-based multi-swarm PSO method for parameter estimation in stochastic biological systems exploiting discrete-time target series. In M. Giacobini, L. Vanneschi, and W. S. Bush, editors, Evolutionary Computation, Machine Learning and Data Mining in Bioinformatics (Proceedings EvoBIO 2012), volume 7246 of LNCS, pages 74–85. Springer, 2012.
50. M. S. Nobile, D. Besozzi, P. Cazzaniga, G. Mauri, and D. Pescini. cupSODA: a CUDA-powered simulator of mass-action kinetics. In V. Malyshkin, editor, Proceedings of 12th International Conference on Parallel Computing Technologies (PaCT 2013), volume 7979 of LNCS, pages 344–357. Springer, 2013.
51. M. S. Nobile, D. Besozzi, P. Cazzaniga, D. Pescini, and G. Mauri. Reverse engineering of kinetic reaction networks by means of cartesian genetic programming and particle swarm optimization. In Evolutionary Computation (CEC), 2013 IEEE Congress on, pages 1594–1601. IEEE, 2013.
52. M. S. Nobile, P. Cazzaniga, D. Besozzi, D. Pescini, and G. Mauri. cuTauLeaping: A GPU-powered tau-leaping stochastic simulator for massive parallel analyses of biological systems. PLoS ONE, 9(3):e91963, 2014.
53. N. Noman and H. Iba. Inference of gene regulatory networks using S-system and differential evolution. In Proceedings of the 2005 Conference on Genetic and Evolutionary Computation, pages 439–446. ACM, 2005.
54. J. Nummela and J. A. Bryant. Evolving Petri nets to represent metabolic pathways. In Proceedings of the 2005 Conference on Genetic and Evolutionary Computation, pages 2133–2139. ACM, 2005.
55. Nvidia. Nvidia CUDA C Programming Guide v5.0, 2012.
56. Nvidia. NVIDIA’s Next Generation CUDA Compute Architecture: Kepler GK110, 2012.
57. J. L. Payne, N. A. Sinnott-Armstrong, and J. H. Moore. Exploiting graphics processing units for computational biology and bioinformatics. Interdisciplinary Sciences, Computational Life Sciences, 2(3):213–220, 2010.
58. B. E. Perrin, L. Ralaivola, A. Mazurie, S. Bottani, J. Mallet, and F. d’Alche Buc. Gene networks inference using dynamic Bayesian networks. Bioinformatics, 19(2):138–148, 2003.
59. L. Petzold. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM Journal of Scientific and Statistical Computing, 4:136–148, 1983.
60. J. W. Pinney, D. R. Westhead, and G. A. McConkey. Petri Net representations in systems biology. Biochemical Society Transactions, 31:1513–1515, 2003.
61. T. Quarles, A. R. Newton, D. O. Pederson, and A. Sangiovanni-Vincentelli. SPICE 3 Version 3F5 User’s Manual. Departement Electrical Engineering Computer Sciences, University of California, 1994.
62. H. Salis and Y. Kaznessis. Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions. The Journal of Chemical Physics, 122(5):054103, 2005.
63. M. Santillán. On the use of the Hill functions in mathematical models of gene regulatory networks. Mathematical Modelling of Natural Phenomena, 3(2):85–97, 2008.
64. M. A. Savageau. Comparison of classical and autogenous systems of regulation in inducible operons. Nature, 252:546–549, 1974.
65. M. A. Savageau. Biochemical Systems Analysis. A Study of Function and Design in Molecular Biology. Addison-Wesley, 1976.
66. R. Storn and K. Price. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 11(4):341–359, 1997.
67. M. Sugimoto, S. Kikuchi, and M. Tomita. Reverse engineering of biochemical equations from time-course data by means of genetic programming. BioSystems, 80:155–164, 2005.
68. G. Szederkenyi, J. R. Banga, and A. A. Alonso. Inference of complex biological networks: distinguishability issues and optimization-based solutions. BMC Systems Biology, 5(177), 2011.
69. M. Vigelius, A. Lane, and B. Meyer. Accelerating reaction-diffusion simulations with general-purpose graphics processing units. Bioinformatics, 27(2):288–290, 2011.
70. J. N. Weiss. The Hill equation revisited: uses and misuses. The FASEB Journal, 11(11): 835–841, 1997.
71. O. Wolkenhauer, M. Ullah, W. Kolch, and C. Kwang-Hyun. Modeling and simulation of intracellular dynamics: choosing an appropriate framework. IEEE Transactions on Nanobiosciences, 3(3):200–7, 2004.
72. S. Xu and Y. Rahmat-Samii. Boundary conditions in particle swarm optimization revisited. IEEE Transactions on Antennas and Propagation, 55(3):760–765, 2007.
73. Y. Zhou, J. Liepe, X. Sheng, M. P. H. Stumpf, and C. Barnes. GPU accelerated biochemical network simulation. Bioinformatics, 27(6):874–876, 2011.