1. A method for establishment of causal regulation network from time series data, comprising the following steps: obtaining a time-series dataset comprising data obtained at a series of time points; generating from a computer tool a first predetermined amount of causal regulation network templates as candidate networks, each candidate network structure being represented by an integer array for genetic algorithm operations; fitting data of said dataset into said candidate networks by downhill simplex method; calculating joint probability density (score) of said fitted candidate networks, S, according to the following equation, respectively: score ( S ) = log ( P ( x  S , Θ ⋒ ) - _ ) ) - d 2 log ( R · n · ( T - 1 ) ) wherein score(S) approximates the logarithm of the joint probability density of a candidate network, S, and the parameters, {circumflex over (Θ)}, represent the values of the parameters Θ that maximize the likelihood function P(x|S,Θ), wherein: P ( x  S , Θ ) = ∏ t = 1 T - 1 ∏ i = 1 n p ( x i , t ; θ i ) , wherein p(x i , t; θ i ) represents conditional probability density of datum i in said candidate network at time t, relatively to parameters θ i , Θ is the cell containing parameters θ l , n represents number of data included in said candidate networks and T represents number of time points, and x i is the value of datum i and is non-negative real; wherein: p ( x i , t + 1 ; θ i ) = N ( x i ( t + 1 ) - α i ∏ j x j ( t ) ω ji - ( 1 - β i ) x i ( t ) , δ 2 ) wherein exponent, ω ji , is real, wherein, when ω ji is positive (negative), datum j induces (represses) the value of datum i and, α i and β i are both positive constants, quantifying the rate of induction and repression, respectively, x i ( t + 1 ) = α i ∏ j x j ( t ) ω ji + ( 1 - β i ) x i ( t ) + ɛ i ( t ) wherein ε i (t) is the measurement error for datum i at time t, δ 2 is variance of distribution in a Gaussian distribution representing said errors, ε i (t); and d 2 log ( R · n · ( T - 1 ) ) _ is a penalty against structure complexity, d equals to the number of parameters in said candidate network, n represents number of data included in said fitted candidate networks, T represents number of time points 0, 1, 2,. .. , T−1, and R represents number of repeats of said measurement; selecting from said fitted candidate networks a second predetermined number of candidate networks with the highest score values; mating and mutating said selected candidate networks according to the genetic algorithm to generate offspring candidate networks; fitting data of said dataset into said offspring networks and calculating score values of said offspring networks; determining whether the average of the obtained score values reaches a plateau; if not, repeating steps of said mating and mutating and said score calculation; otherwise, selecting a third predetermined number of the offspring candidate networks with highest score values; and outputting said selected offspring candidate networks as result of establishment of causal regulation network.