From 2009.igem.org

Genetic Algorithm (GA)

Introduction to GA

What's GA

Genetic algorithms are powerful methods for complex optimization problems. They are essentially evolution in a computer. A population of GA population random candidate solutions is generated and scored. These solutions are allowed to mate with each other, generating offspring that are hybrids of both parental solutions. The mating frequencies are related to a boltzmann-like probability term; akin to simulated annealing, the simulation temperature required for calculating these is high early in the run, and cools down. Since the mating probability for a given candidate solution is related to its score, better solutions mate more frequently, making their components more numerous in the offspring. This new population of candidate solutions is scored, ranked, and allowed to mate.

Why GA?

Genetic algorithm (GA) is utilized to identify the topology of the network. A rate term can be described as two parts: signal and function type. The former usually indicates the positive and negative interaction from reactants to products which strongly correlate with the SBGN entity relationship diagram. The later reflects the biochemical reaction type. GA searches an optimized solution by mimicking the process of evolution. Members of an initial population compete for the privilege to mate and produce offspring for successive generations. The probability of an individual mating is proportional to the fitness of the environment. The generation with better fitness will be obtained generation by generation. As the topology, the interaction forms and the parameters for each component function are all necessary in the final result, we employed a multi-level optimization. The lowest level optimization is to give the suitable parameters with the topology and the interaction forms known. For the next level, we are trying to find the interactions forms of component functions with the known topology with genetic algorithm. At last, for the upper level, the potential topologies are selected also with genetic algorithm. And the fitness functions are the similarity between the dynamic behaviors and the target, which are measured by the standard error.

Optimization Processes of GA

Step1: Randomly generate a population of topologies.

Step2: For each topology, generate a population of networks with a serial of random interaction forms.

Step3: For each network, obtain the fitness score by optimize the parameters.

Step4: Obtain the fitness score of a topology by searching the most favorable network.

Step5: Evolutes the topologies with each fitness score.

Parameter For GA

Population: the number of members in one generation

Mutation ratio: the ratio of each gene alter its status in a generation

Recombine ratio: the ratio that whether a gene change with its allel in mating

Max_cycle: the maximum number of generations considered for the evolution

Particle Swarm Optimization Algorithm (PSO)

Introduction to PSO

What's PSO

Particle Swarm Optimization(PSO) is an optimization algorithm based on swarm intelligence theory. Motivated by the evolution of nature, a series of evolutionary computation techniques, such as evolutionary programming, genetic algorithms, evolutionary strategies, are proposed, in which a best solution is evolved through the generations. In contrast to evolutionary computation techniques, Eberhart and Kennedy developed a different algorithm through simulating social behavior of bird flocking or fish schooling. In a particle swarm optimizer, instead of using genetic operators, individual swarms are evolved through cooperation and competition. Each particle adjusts its flying according to its own flying experience and its companions' flying experience. The position of each particle represents a potential solution to a problem.

The following is excerpt from"http://www.swarmintelligence.org/", it shows clearly how PSO works: "Each particle keeps track of its coordinates in the problem space which are associated with the best solution (fitness) it has achieved so far. (The fitness value is also stored.) This value is called pbest. Another "best" value that is tracked by the particle swarm optimizer is the best value, obtained so far by any particle in the neighbors of the particle. This location is called lbest. when a particle takes all the population as its topological neighbors, the best value is a global best and is called gbest.

The particle swarm optimization concept consists of, at each time step, changing the velocity of (accelerating) each particle toward its pbest and lbest locations (local version of PSO). Acceleration is weighted by a random term, with separate random numbers being generated for acceleration toward pbest and lbest locations." (cite from the website)

Why PSO?

In past several years, PSO has been successfully applied in many research and application areas. It is demonstrated that PSO gets better results in a faster, cheaper way compared with other methods. Compared to GA, the advantages of PSO are that PSO is easy to implement and there are few parameters to adjust. The most important reason we choose to implement this algorithm into our software is that this algorithm is easy to realize parallelization. Since the most time-consuming part in our scheme is the optimization of parameters for a given topological structure. If we cannot find a efficient optimizer, it is impossible to deal with systems contains more than five or six nodes. parallelization of the optimization process will be implemented in our next version.

Basic Formulas of PSO

First we must be clear with several concepts:

particle:have a position and a velocity, fly through parameter space

velocity: velocity of a particle in the parameter space

pbest: The best solution a particle has achieved so far

lbest: Best solution obtained so far by any particle in the neighbors of the particle

gbest: The best solution achieved by all particles so far

Parameters used in PSO algorithm:

Number of Particles:

learning factors <math>c_{1}</math>, <math>c_{2}</math>

inertial factor <math>w</math>

range of parameter and velocity on each dimension

optimization threshold

The particles update their positions and velocities with following equation:

<math>v[] = w*v[] + c_{1}*rand()*(pbest[]-postion[]) + c_{2}*rand()*(gbest[] - position[])</math>

<math>position[] = position[] + v[]</math>

where <math>c_{1}</math> and <math>c_{2}</math> are two learning factors, we choose <math>c_{1}=c_{2}=2</math>. Particles' velocities on each dimension are clamped to a maximum velocity <math>V_{\max }</math>. If the sum of accelerations would cause the velocity on that dimension to exceed <math>V_{\max }</math>. Then the velocity on that dimension is limited to <math>V_{\max }</math>. If a particle's position go beyond the range of parameter, we will set the parameter on that dimension to the boundary condition and let the sign of velocity on that dimension change meanwhile in order to simulate reflection.

Performance of PSO

In contrast with simulated annealing, another method we used to optimize parameters, PSO performs much better both on the efficiency and accuracy. Specifically, PSO is good at dealing with very large parameters because it is a global optimization algorithm.

Reference

http://www.swarmintelligence.org/

Eberhart, R.C. and Kennedy, J (1995). A new optimizer Using Particle Swarm Theory, Proc. Sixth International Symposium on Micro Machine and Human Science

Data Structure and Organization

struct reactor{
	int rnum;//number of reactants
	char **ract;//names of reactants
};

struct product{
	int pnum;//number of products
	char **prdt;//names of products
};

struct reaction{
	int num;//reaction number
	struct reactor RACT;//info of reactants
	struct product PRDT;//info of products
	int type;//reaction type as listed in keneticlaw.pdf
	double *para;//all parameters
	char note[MAXLEN];//self-defined notes
};

Global Sensitivity Analysis (GSA)

A general objective in our modeling process is an exploration of the high-dimensional input variable space as thoroughly as possible for its impact on observable system behavior, often with either optimization in mind or simply for achieving a better understanding of the phenomena involved. Since the system input->output behavior is typically a nonlinear relationship, simple logic suggests that the number of runs could grow exponentially with the number of input variables. However, an emerging family of high dimensional model representation concepts and techniques capable of dealing with such input->output problems in a practical fashion. RS-HDMR is a global sensitivity analysis technique that can decompose the high-dimensional, nonlinear contributions of each parameter to the network properties(represented by their total sensitivity) into a hierarchy of low-dimensional terms. Calculations by RS-HDMR require only the ODE model, an estimate of the initial conditions and an estimate of the dynamic range to explore for each parameter.

Introduction to GSA

What's GSA?

Why GSA?

Theoretical Foundations of GSA

High dimensional model representation(HDMR) is a general set of quantitative model assessment and analysis tools for capturing high dimensional IO system behavior. As the impact of the multiple input variables on the output can be independent and cooperative, it is natural to express the model output <math>f(x)</math> as a finite hierarchical correlated function expansion in terms of the input variables:

Where the zeroth-order component function <math>f_{0}</math> is a constant representing the mean response to <math>f(x)</math> , and the first order component function <math>f_{i}(x_{i})</math> gives the independent contribution to <math>f(x)</math> by the ith input variable acting alone, the second order component function <math>f_{ij}(x_{i},x_{j})</math> gives the pair correlated contribution to <math>f(x)</math> by the input variables <math>x_{i}</math> and <math>x_{j}</math>, etc. The last term contains any residual nth order correlated contribution of all input variables. The basic conjecture underlying HDMR is that the component functions arising in typical real problems are likely to exhibit only low order l cooperativity among the input variables such that the significant terms in the HDMR expansion are expected to satisfy the relation: <math>l<<n</math> for <math>n>>l</math> . An HDMR expansion to second order

often provides a satisfactory description of <math>f(x)</math> for many high dimensional systems when the input variables are properly chosen. This is also the formula we use to achieve HDMR expansion. How to determine HDMR component functions? Practical formulations of the HDMR component functions determine whether we can carry out this algorithm. Theoretically, a component function

is obtained by minimizing the functional,

<math>\min_{f_{i_{1}i_{2}...i_{l}}(u_{i1},u_{i2},...,u_{il}}\int_{\Omega}\omega_{i_{1}i_{2}...i_{l}}(\hat{x},u) [f(u)-f_{0}-\sum_{i=1}^{n}f_{i}(u_{i})-\sum_{1\leq i<j\leq n} f_{ij}(u_{i}, u_{j})-...- \sum_{1\leq i_{1}<...<i_{l}\leq n}f_{i_{1}i_{2}...i_{l}}(u_{i1},u_{i2},...,u_{il})]^{2}du</math>

under a suitable specified orhogonality condition which guarantees that all the component functions are determined step by step. Here,<math>\hat{x}=(x_{i_{1}},x_{i_{2}},...,x_{i_{n}}</math> , <math>\omega_{i_{1}i_{2}...i_{l}}(\hat{x},u)</math> is a weight function. Different weight functions will produce distinct, but formally equivalent HDMR expansions. There are two commonly used HDMR expansions: Cut- and RS(Random Sampling)-HDMR. Cut-HDMR expresses <math>f(x)</math> in reference to a specified cut point, while RS-HDMR depends on the average value of <math>f(x)</math> over the whole domain <math>\Omega</math> . We adopt RS-HDMR to do global sensitivity analysis. For RS-HDMR, we first rescale variables <math>x_{i}</math> such that <math>0<x_{i}<1</math> for all i. The output function <math>f(x)</math> is then defined in the unit hypercube

by suitable transformations. The component functions of RS-HDMR possess the following forms:

Where <math>du^{i}</math> and <math>du^{ij}</math> are just the product <math>du_{1}du_{2}...du_{n}</math> without <math>du_{i}</math> and <math>du_{ij}</math> respectively. Evaluation of the high dimensional integrals in the RS-HDMR expansion is carried out by the Monte Carlo random sampling. Orthogonality conditions are used to obtain the above formulas. The integral of a component function of RS-HDMR with respect to any of its own variables is zero, i.e.,

Using the orthobonality property of the RS-HDMR component functions, it can be proven that the total variance of <math>f(x)</math> caused by all input variables sampled uniformly over their full range may be decomposed in to distinct input contributions in the following manner.

<math>\begin{align} \sigma_{\bar{f}}^{2} &= \int_{K^{n}}[f(x)-\bar{f}]^{2}dx \\ &= \sum_{i=1}^{n}\int_{0}^{1}f_{i}^{2}(x_{i})dx_{i}+\sum_{1\leq i< j \leq n} \int_{0}^{1}\int_{0}^{1}f_{ij}^{2}(x_{i},x_{j})dx_{i}dx_{j}+... \\ &= \sum_{i=1}^{n}\sigma_{i}^{2}+\sum_{1\leq i< j \leq n}\sigma_{ij}^{2}+... \\ \end{align}\,\!</math>

Thus, the total variance <math>\sigma_{\bar{f}}^{2}</math> is the sum of first-order variances <math>\sigma_{i}^{2}</math> , second-order covariances <math>\sigma_{ij}^{2}</math> , etc. The magnitudes of the indices <math>\sigma_{i}^{2}</math>, <math>\sigma_{ij}^{2}</math> reveal how the output uncertainty is influenced by the input uncertainties and the nature of the cooperativities that exist. The direct determination of the component functions of RS-HDMR at different values in the unit hypercube by Monte Carlo integration can require a large number of random samples. To reduce the sampling effort, the RS-HDMR component functions may be approximated to any desired level of accuracy by either analytical approximation or numerical approximation. In our code, we implement the former approximation, which approximate the RS-HDMR component function by expansion in terms of a suitable set of functions. Specifically, orthogonal polynomials \cite{paper7} are used as a basis to approximate <math>f_{i}(x_{i})</math> , <math>f_{ij}(x_{i},x_{j})</math> as follows:

<math>f_{i}(x_{i})\sum_{k=1}^{\infty}\alpha_{k}^{i}\phi_{k}(x_{i})</math>

<math>f_{ij}(x_{i},x_{j})=\sum_{k,l=1}^{\infty}\beta_{kl}^{ij}\phi_{k}(x_{i})\phi_{l}(x_{j})</math>

In most cases, to achieve a desired accuracy using <math>\phi_{1}(x)</math> ,<math>\phi_{2}(x)</math> and <math>\phi_{3}(x)</math> is sufficient. Description of the use of orthonormal polynomials in details can be found in ref\cite{paper7}, we only give some results here.

<math>\begin{align} \sigma_{i}^{2} &= \int_{0}^{1}f_{i}^{2}(x_{i})dx_{i}\approx \int_{0}^{1}[\sum_{k=1}^{s_{1}}\alpha_{k}^{i}\phi_{k}(x_{i})]dx_{i} \\ &= \sum_{k=1}^{s_{l}}(\alpha_{k}^{i})^{2} \\ \end{align}\,\!</math>

<math>\begin{align} \sigma_{ij}^{2} &= \int_{0}^{1}\int_{0}^{1}f_{ij}^{2}(x_{i},x_{j})dx_{i}dx_{j} \\ &\approx \int_{0}^{1}\int_{0}^{1}[\sum_{k=1}^{s_{1}}\sum_{k=1}^{s_{1}'}\beta_{kl}^{ij}\phi_{k}(x_{i})\phi_{l}(x_{j})]^{2}dx_{i}dx_{j} \\ &= \sum_{k=1}^{s_{1}}\sum_{k=1}^{s_{1}'}(\beta_{kl}^{ij})^{2}\\ \end{align}\,\!</math>

Only one random sampling set of <math>f(x)</math> is needed to determine <math>f_{0}</math> and all the expansion coefficients and consequently,<math>\sigma_{\bar{f}}^{2}</math>, <math>\sigma_{i}^{2}</math> , <math>\sigma_{ij}^{2}</math> . This dramatically reduces the sampling effort and makes global sensitivity analysis of a model very efficient.

Local Sensitivity Analysis (LSA)

Besides global sensitivity analysis, we also implement additional local sensitivity analysis technique, which can measure time-varying sensitivities along arbitrary trajectories. For local methods, the basic measures of sensitivity are the partial derivatives.

<math>\frac{\partial(Model Responses)}{\partial (Model Parameters)}</math>,

called the sensitivity coefficients of the model\cite{paper8}. We adopt the method proposed by \cite{paper9}, which can determine the sensitivity to the perturbation throughout the time evolution of the system, regardless of the nature of the trajectory. The network is modeled by the ordinary differential equation

<math>\frac{d}{dt}\mathbf{s}(t)=N\mathbf{v}(\mathbf{s}(t),\mathbf{p},t)</math>

where the n-vector s is composed of the concentrations of each species, the constant r-vector p is composed of the parameters of interest in the model, the m-vector valued function

<math>\mathbf{v}=\mathbf{v}(s, p, t)</math>

describes the rate of each reaction as function of species concentrations and parameter values.

Introduction to LSA

What's LSA?

Why LSA?

Basic Definitions

Definition 1 Given an initial condition <math>s(0)=s_{0}</math> and a set of parameter values <math>p_{0}</math> , which together for a vector <math>q_{0}</math> ,the time-varying concentration sensitivity coefficients(or concentration response coefficients) are define as the elements of the <math>n\times(n+r)</math> matrix function <math>R_{q}^{s}(t)</math> given by

<math>R_{q}^{s}(t):=\frac{\partial s(t,q)}{\partial q}|_{q=q_{0}}=\lim_{\Delta q\rightarrow0} \frac{s(t,q_{0}+\Delta q)-s(t,q_{0})}{\Delta q}</math>

for all <math>t\geq0</math> The response coefficient defined above provides a measure of the difference between this "perturbed trajectory" and the "nominal" trajectory at each time t. Definition 2 Given an initial condition <math>s(0)=s_{0}</math> and a set of parameter values <math>p_{0}</math> , which together for a vector <math>q_{0}</math> ,the time-varying rate sensitivity coefficients(or rate response coefficients) as the elements of the <math>m\times r</math> matrix <math>R_{q}^{v}(t)</math> function given by,

<math>R_{q}^{v}(t)=\frac{\partial v(s(t,q),p,t)}{\partial q}|_{q=q_{0}}</math>

for all, where the derivatives are evaluated at and The rate response coefficients give the response in the rates at time t to a perturbation at time zero.

Computation

The sensitivity coefficients are defined by a first-order linear ordinary differential equation,

<math>\frac{d}{dt}\frac{\partial s(t,q)}{\partial q}=N[\frac{\partial v(t)}{s}\frac{\partial s(t)}{\partial q}+ \frac{\partial v(t)}{\partial q}]</math>

for all <math>t \geq 0</math>. The choice of <math>s(0)=s_{0}</math> and <math>p_{0}</math> fix a trajectory <math>s(t)=s(t,s_{0},p_{0}</math> , and hence determines the function <math>v(t)=v(s(t),p,t)</math> as well. The equation above can be solved numerically for the concentration sensitivities. In performing such a calculation, we must note that the right-hand side depends on the time-varying value of the species concentration vector s(t), so a simple computational strategy is to determine <math>s_{i}(.)</math> and <math>\partial s_{i}(.)/\partial q</math> simultaneously.

Team:USTC Software/Algorithm