Technical Description of CADLIVE SIMULATOR
A. PARSEDAE............................................................................................................................................................. 6
1 Abstract ................................................................................................................................................................. 6
1.1 Input data........................................................................................................................................................ 6
1.2 Mathematical formulas................................................................................................................................... 7
2 Input data............................................................................................................................................................... 7
2.1 Structure of input data .................................................................................................................................... 7
2.2 Species naming rules...................................................................................................................................... 9
2.3 Complex ......................................................................................................................................................... 9
2.4 Species Attribute .......................................................................................................................................... 10
2.5 Definition of regulator-reaction equations ................................................................................................... 10
3 Conversion to mathematical models ................................................................................................................... 13
3.1 TT (ordinary Transcription and Translation) ................................................................................................ 13
3.1.1 Transcription without regulators (gene0) .............................................................................................. 14
3.1.2 Transcription with regulators (gene1) ................................................................................................... 14
3.1.3 Translation............................................................................................................................................. 14
3.2 CMA (Conventional Mass Action)............................................................................................................... 14
3.2.1 Binding reactions................................................................................................................................... 14
3.2.2 Conversion reaction............................................................................................................................... 15
3.2.3 General enzyme reaction ....................................................................................................................... 15
3.2.4 Transport ............................................................................................................................................... 16
3.2.5 Spontaneous decomposition .................................................................................................................. 16
3.3 GMA (General Mass Action) ....................................................................................................................... 17
3.4 MM (simplified Michaelis-Menten equations) ............................................................................................ 18
3.5 DAEs ............................................................................................................................................................ 18
3.5.1 Steady-state approximation I................................................................................................................. 18
3.5.2 Steady-state approximation II................................................................................................................ 19
3.5.3 Rapid equilibrium approximation ......................................................................................................... 19
3.5.4 Complexes in TPP ................................................................................................................................. 20
3.6 S-system ....................................................................................................................................................... 20
3.7 Connection between layers........................................................................................................................... 21
4. Selection of conversion methods........................................................................................................................ 22
4.1 Gene-protein layer........................................................................................................................................ 22
4.2 Metabolic layer............................................................................................................................................. 22
4.3 Connection between the layers..................................................................................................................... 22
B. CHECKDAE ......................................................................................................................................................... 23
1
KIT, Kurata Laboratory

1 Objectives............................................................................................................................................................ 23
2. Function.............................................................................................................................................................. 23
2.1 Grammar check for mathematical models.................................................................................................... 23
2.2 Expansion to mathematical equation trees ................................................................................................... 24
2.3 Function for checking parameters and variables .......................................................................................... 24
2.4 Function for expanding temporary expressions ........................................................................................... 24
2.5 Function for indexing variables and parameters .......................................................................................... 24
2.6 Function for replacing the indexes ............................................................................................................... 24
2.7 Symbolic partial differentiation.................................................................................................................... 24
2.8 Function for S-system conversion................................................................................................................ 24
2.9 Parameter survey .......................................................................................................................................... 24
3 Specification of input and output ........................................................................................................................ 24
3.1 Specification of input ................................................................................................................................... 24
3.2 Specification of output ................................................................................................................................. 24
3.2.1 User function file................................................................................................................................... 24
3.2.2 Parameter file ........................................................................................................................................ 25
3.2.3 Temporary mathematical expression file............................................................................................... 27
3.2.4 S-system coefficient file ........................................................................................................................ 27
4 Data structure ...................................................................................................................................................... 28
4.1 Main structures............................................................................................................................................. 28
4.2 Node_T structure.......................................................................................................................................... 29
5 Procedure............................................................................................................................................................. 30
5.1 Input procedure (Read_data) ........................................................................................................................ 33
5.1.1 Reading checkdae files: data split and data number acquisition (split_dae_file) .................................. 33
5.1.2 Index check and memory allocation (allocate_area). ............................................................................ 33
5.1.3 Reading parameter labels (read_param_base)....................................................................................... 33
5.1.4 Reading parameters (read_param)......................................................................................................... 33
5.1.5 Reading variables (read_var)................................................................................................................. 33
5.1.6 Reading temporary mathematical expressions (read_mid) ................................................................... 33
5.1.7 Reading differential equations and algebraic equations (read_equation) .............................................. 33
5.2 Expansion of mathematical equations in a tree structure (set_expression、Make_node) ........................... 34
5.2.1 Reading mathematical expression (get_token)...................................................................................... 34
5.2.2 Checking mathematical expression (check_token) ............................................................................... 34
5.2.3 Adding a blank child node (add_emptyChild)....................................................................................... 34
5.2.4 Inserting a child node (insert_childNode) ............................................................................................. 34
5.2.5 Inserting a parent node (insert_parentNode) ......................................................................................... 34
5.2.6 Deleting nodes (del_node)..................................................................................................................... 34
5.2.7 Searching equivalent or less prior nodes (seek_same) .......................................................................... 35
2
KIT, Kurata Laboratory5.2.8 Reconstructing mathematical expressions (Set_expr)........................................................................... 35
5.3 Checking variables (check_node、check_detail) ........................................................................................ 35
5.4 Expanding temporary mathematical expressions (expand_mid).................................................................. 35
5.5 Assigning an index to a variable (set_variable_index)................................................................................. 35
5.6 Converting variables (set_numerical_index, check_node, replace_variable) .............................................. 35
5.7 Symbolic partial differentiation (Make_diff) ............................................................................................... 35
5.8 Output process.............................................................................................................................................. 36
5.8.1 Output of user functions for simulators (Write_userfunc)..................................................................... 36
5.8.2 Output of parameter files (Write_paramFile) ........................................................................................ 36
5.8.3 Output of temporary mathematical expression obtaining file

(Write_flux) ....................................... 36

5.9 S-system conversion（Set_Ssystem） ........................................................................................................ 36
3.2.5 S-system user function file .................................................................................................................... 37
C. SOLVER................................................................................................................................................................ 39
1. Introduction ........................................................................................................................................................ 39
2. Analysis types..................................................................................................................................................... 39
3. Input for control data.......................................................................................................................................... 39
D. MERGEPARAM ................................................................................................................................................... 41
1. Objectives........................................................................................................................................................... 41
2. Function.............................................................................................................................................................. 41
3. Input/output specification................................................................................................................................... 41
3.1 Y_START ..................................................................................................................................................... 41
3.2 PARAM ........................................................................................................................................................ 41
4. Command ........................................................................................................................................................... 42
5. Error message..................................................................................................................................................... 42
E. SENSITIVITY AND STABILITY ANALYSIS BY S-SYSTEM .......................................................................... 43
1. Objectives........................................................................................................................................................... 43
2. Introduction of S-system .................................................................................................................................... 43
2.1 S-system conversion..................................................................................................................................... 43
2.2 Sensitivity analysis....................................................................................................................................... 44
2.2.1. Logarithmic gain of a metabolite ......................................................................................................... 44
2.2.2. Logarithmic gain of a flux.................................................................................................................... 45
2.2.3. Sensitivity with respect to a rate constant ............................................................................................ 45
2.2.4. Sensitivity of dependent variables with respect to a kinetic order ....................................................... 45
2.2.5. Sensitivity of a flux to a rate constant .................................................................................................. 46
2.2.6. Sensitivity of a flux with respect to a kinetic order.............................................................................. 46
2.3 Stability analysis .......................................................................................................................................... 46
3 Process flow ........................................................................................................................................................ 47
4 Function............................................................................................................................................................... 49
3
KIT, Kurata Laboratory4.1 Creation of S-system parameter file ............................................................................................................. 49
4.2 Sensitivity analysis....................................................................................................................................... 49
4.3 Stability analysis .......................................................................................................................................... 49
5 Input/Output Specification .................................................................................................................................. 49
5.1 Input specification ........................................................................................................................................ 49
5.2 Output specification ..................................................................................................................................... 50
F. OPTIMIZER (GA) ................................................................................................................................................. 52
1. Introduction ........................................................................................................................................................ 52
2. Function specification and application problems ............................................................................................... 52
2.1 Function Specification.................................................................................................................................. 52
2.2 Application problem..................................................................................................................................... 53
3. Execution of programs ....................................................................................................................................... 53
4. Input/output specification................................................................................................................................... 54
4.1 Parameter-range setting file.......................................................................................................................... 54
4.2 GA-parameter setting file............................................................................................................................. 55
4.2.1 GA-parameter setting file...................................................................................................................... 55
4.2.2 Examples for GA parameter setting file ................................................................................................ 57
4.2.3 GA-population setting file..................................................................................................................... 58
4.2.4 Examples for GA-population setting file .............................................................................................. 58
4.3 Standard output ............................................................................................................................................ 59
4.4 Output file................................................................................................................................................... 61
5. Detailed specification of functions..................................................................................................................... 61
5.1 Encode method............................................................................................................................................. 61
5.1.1 Bit-string GA......................................................................................................................................... 61
5.1.2 Real GA................................................................................................................................................. 62
5.2 GA type ........................................................................................................................................................ 62
5.2.1 Distributed GA ...................................................................................................................................... 62
5.2.2 Distributed and integrated GA (DIGA) ................................................................................................. 62
5.3 Generation alternation .................................................................................................................................. 62
5.3.1 Ordinary generation alternation method................................................................................................ 63
5.3.2 Minimal Generation Gap (MGG).......................................................................................................... 63
5.4 Crossover method......................................................................................................................................... 63
5.4.1 BLX-α [RGA] ....................................................................................................................................... 63
5.4.2 UNDX / UNDX-m [RGA] .................................................................................................................... 64
5.4.3 SPX [RGA] ........................................................................................................................................... 65
5.4.4 N point crossover [BGA] ...................................................................................................................... 66
5.5 Mutation method .......................................................................................................................................... 66
5.5.1 Uniform mutation within the region [RGA].......................................................................................... 66
4
KIT, Kurata Laboratory5.5.2 Uniform mutation with fixed width [RGA]........................................................................................... 66
5.5.3 Normal mutation with fixed width [RGA] ............................................................................................ 66
5.5.4 Uniform mutation with variable width [RGA] ...................................................................................... 67
5.5.5 Normal mutation with variable width [RGA] ....................................................................................... 67
5.5.6 Bit reverse mutation [BGA] .................................................................................................................. 67
6. Flow chart........................................................................................................................................................... 67
6.1 Single without MPI..................................................................................................................................... 67
6.2 MPI............................................................................................................................................................... 68
6.2.1 Master-slave model ............................................................................................................................... 68
6.2.2 A flow char in MPI ................................................................................................................................ 69
G. CLIENT-SERVER MODEL .................................................................................................................................. 72

5
KIT, Kurata Laboratory

A. PARSEDAE
1 Abstract
In order to convert the regulator-reaction equations into mathematical models, CADLIVE classifies biochemicalractions into three layers: gene, protein, and metabolic layers, and divides the conversion process into two stages,
the first conversion (ordinary transcription and translation equations = TT, conventional mass action = CMA,
general mass action = GMA, simplified Michaelis-Menten equations = MM), and the second conversion(differential and algebraic equations = DAEs, S-system). From the standpoint of mathematical conversions, theapplied mathematical conversion strongly depends on the layer that the regulator-reaction equations belong to. At
the first stage, both gene and protein networks employ TT and CMA, whereas the metabolic network uses GMA or
MM. In the gene layer, since various molecules such as proteins, amino acids, nucleic acids, and RNAs act inconcert for transcription and translation, it is difficult to mathematically describe such reactions based on theirconcrete molecular mechanism. The use of TT is a rational choice for taking in gene regulations within a cell. In theconversion of metabolic networks into GMA or MM, the concentrations of enzyme-metabolite complexes arecancelled compared with those of metabolites, because the former are far less than the latter. By contrast, in theprotein layer, many proteins function in a complex or modified form, thus it is not practical to cancel theconcentrations of the active complexes or modified ones. Thus, the use of CMA and TT describes the protein signaltransduction pathways. Problems for CMA are that its differential equations are stiff due to the huge differences invalues of kinetic parameters and molecular concentrations, and that many kinetic parameters are required.
At the second stage, in order to overcome such problems, we applied the Two-Phase Partition (TPP) method tothe conversion of CMA with TT to differential and algebraic equations (DAEs), thereby reducing not only thestiffness, but also decreasing the number of kinetic parameters. The TPP method divides the kinetics of molecularinteractions into two phases, the molecular binding phase and the reaction phase, assuming thatassociation/dissociation rates between proteins to be quite fast compared with the rates of synthesis/degradation ofmRNAs or proteins. This is a commonly used assumption in biological systems. Consequently, TPP substitutes
algebraic equations for stiff differential equations. In addition to DAEs, CADLIVE is able to convert the ordinary
differential equations of TT, CMA, GMA, and MM into S-system at the steady state, which is useful for analyzing
the sensitivity and stability in symbolic form.
The parsedae module converts the regulator-reaction equations (sanac file Fig. A1), which are constructed by the
CADLIVE editors, into the mathematical equations (checkdae file) including TT, CMA, GMA, MM, DAEs, andS-system. The checkdae file employs the indexes of species' names so that one can understand the model or edit it
manually. One is allowed to edit the checkdae file directly according to the instruction. The checkdae file is written
in the text format.

1.1 Input data
The parsedae module inputs the "sanac" file (Fig. A1) that mainly consists of two kinds of data. One is the species(molecules) with their associated attributes. The other is the regulator-reaction equation that determines theinteraction among species.
6
KIT, Kurata Laboratory1.2 Mathematical formulas
The following mathematical formulas have been employed.
Protein layer
CMA (Conventional Mass Action) is applied to the protein layer, focusing on the formation of complexes andmodified proteins.
Gene layer
TT (ordinary transcription-translation equations) is applied to the gene layer.
Metabolic layer
GMA (General Mass Action) is applied to the metabolic layer, neglecting detailed mechanisms including theformation of enzyme-metabolite complexes.
MM (simplified Michaelis-Menten equations) is applied to the metabolic layer, canceling the concentrations ofenzyme-metabolite complexes. MM is obtained by simplifying ordinary Michaelis-Menten equations.
All the layers
DAEs (differential and algebraic equations) are applied to the gene-protein layer using the two-phase partition(TPP) method.
S-system can be applied to all the layers.

2 Input data
2.1 Structure of input data
Systems Biology Markup Language (SBML) is one of the most advanced markup language that describes
biochemical networks at a concrete molecular interaction level. We extend SBML level 2 to establish the sanac
format that includes all the information necessary for mathematical conversion and dynamic simulation. The
original sanac file (Fig. 1A) contains various data including the coordinates of the species, which are used fordrawing a biochemical network map by the CADLIVE GUI editors. Here we omit the data that are not necessary
for mathematical conversion and dynamic simulation.

7
KIT, Kurata Laboratorysanac1 model (name, modelInfo)
1 listOfCompartments (numberOfCompartments)+ compartment (name, volume)
1 listOfSpecies (numberOfSpecies)+ Specie (specieID, name, compartment)
1 specieAttribute (specieClass, massBalance, totalAmount,decomposition, initialAmount)? listOfElements (numberOfElements):for a complex+ specieReference (specieID, name, compartment, stoichiometry)+ listOfReactions(layer, numberOfReactions)+ Reaction (reactionID, groupID, reactionType, isConvertible,
modifierReactionEquation)
1 listOfReactants (numberOfReactants)* specieReference (specieID, name, compartment, stoichiometry)? listOfProducts (numberOfProducts):for degradation, an empty is possible.

+ specieReference (specieID, name, compartment, stoichiometry)? listOfModifiers (numberOfModifiers)
1 modifieSpecieReference (specieID, name, compartment, stoichiometry,modifierType, isTotal)? modifierComplex
1 specieReference (specieID, name, compartment, stoichiometry)Fig. A1 Conversion rules implemented in the "sanac" representation. They are necessary and sufficient forcomputationally converting biochemical networks into mathematical models. The regulator is exactly the sameas the SBML-defined modifier. The net regions are the extension from SBML level 2. Legend: 1: an elementthat appears once. *: an element that repeats more than zero. +: an element that repeats more than one. ?: anelement that appears once or zero. t: an element that has text in content. (): provided as attributes. The net
regions are the extension from SBML level 2.
The major extensions are as follows:
The attribute of "specieID" refers from the lists of reaction equations or species.
The element of <listOfElements> shows the components of the complexes, which are defined as"specieReference".
The attribute of "layer" is added to the element of <listOfReactions>, whereby a biochemical network canbe divided into three layers, metabolic (metabolic), protein signal transduction (protein), and gene
regulatory networks (gene).
The attribute of "modifierReactionEquation" is added to the element of <reaction>.
8
KIT, Kurata Laboratory

The elements of <listOfModifiers> and <modifierComplex> are added to the element of <Reaction>, which
are necessary for mathematical conversion. The modifierComplex corresponds to substrate-enzyme complex
in ordinary enzyme reactions.
The attribute of "reactionType" is added to distinguish regulator-reaction equations, which are indispensable
for mathematical conversion.
The attribute of "modelInfo" is added to comment the model.

2.2 Species naming rules
The following characters are allowed to name species
Initial character- upper and lower alphabets[a-z] | [A-Z]- number[0-9]- left parentheses(

Second characters- upper and lower alphabets [a-z] | [A-Z]- underbar_- number[0-9]- left and right parentheses- colon:- hyphen_- period.

()

The period may be used for indicating a decimal point.

2.3 Complex
A complex is defined by listing its components as specieReference in <listOfElements>. The use of "stoichiometry"of "specieReference" is able to determine the ratio of components. Not only a monomer/modified but also acomplex is allowed to list as the components of the complex, thus users are able to define a higher complex byusing complex, modified, or monomer. It is required to keep the naming rules. The employed species are required
not to show reference circulation and to be defined in the model.
The element of <listOfElements> is used:
1: to assign the species decomposition to the respective differential equations,
2: to generate the mass balance equations and the differential equations for the total amounts of the species whoseclass is protein or modified in the TPP conversion,
3: to calculate the total concentrations of promoters and enhancers in the "gene1" reaction.
9
KIT, Kurata Laboratory2.4 Species Attribute
The species is defined by three indispensable attributes:
specieID:The integer unique to the combination of species' name with its compartment

name:Users can name the species arbitrarily.

compartment:Location of the species

The other attributes are provided as shown in Table A2.
Table A2 Attributes of species (molecules)
Attribute
specieID

Name

Content
Value
the integer unique to the combinationofspecies'
name
with
its
compartment
name of species
named by user

Compartment
SpecieClass

location of species
kinds of species

massBalance

total mass balance equationis required.

totalAmount

total concentration

decomposition

spontaneous degradation in vivo

initialAmountinitial concentration of species

see text
DNA_gene
DNA_promoter
DNA_enhancer
DNA_othersRNA
proteinmetabolite
environmental_factor
ion_signal
complexmodifier_complex
modified
small_molecule
text_option
otherson
offconstant
variableonoffreal value

2.5 Definition of regulator-reaction equations
A model has multiple elements of <listOfReactions>, which has two indispensable attributes of "layer". The
attribute of "layer" has three values, "metabolic", "protein", and "gene", which indicate where reactions occur. The
metabolic layer indicates metabolic networks, the protein layer signal transduction pathways, and the gene layertranscription and translation. The mathematical formulas are determined by the layer that reaction belongs to.
10
KIT, Kurata Laboratorylayer:determining the layer that reactions occur, which consists of "metabolic", "protein",and "gene".
numberOfReactions:the number of reactions.

The element of <reaction> consists of the elements of <listOfReactants>, <listOfProducts>, <listOfModifiers>,and <modifierComplex>
listOfReactants:defines reactant(s)listOfProducts:defines product(s), For "reactionType = degradation", an empty is possible.

listOfModifiers:define modifier.

modifierComplex:define the complex containing a modifier.

The " reactionType" attribute of <reaction> is shown in Table A3.
Table A3 Values of the attribute of reactionType
Valuebindingbinding_with_stoichiometric_changes
homo_association_or_modification
homo_association_or_modification_with_stoichiometric_changes
elimination
elimination_with_stoichiometric_changes
reversible_conversion
irreversible_conversion
reversible_conversion_regarding_multicomponent
irreversible_conversion_regarding_multicomponent
transport
option_transporttranscription
translation
degradation
The elements of <listOfReactants>, <listOfProducts>, <modifierComplex>, and <listOfElements> have to bereferred from <specieReferences>, which consists of four attributes, "specieID", "name", "compartment", and"stoichiometry".
specieID

:

name

:

compartment:cytoplasm, nucleoplasm, membrane, environment, …..

stoichiometry:defined by real values. The default value is one.

11
KIT, Kurata Laboratory

The element of <listOfModifiers> must be referred from <modifierSpecieReference>, which is the element of<specieReference> with the attributes of "modifierType" and "isTotal".
modifierType:enzyme, activator, inhibitor

isTotal:false, true

Here we explain the detailed rules by using the nitrogen assimilation system. The attribute of "isConvertible" inthe element of <reaction> determines if the regulator-reaction equations can be converted into mathematicalequations. The parsedae module of CADLIVE converts all the reactions with "isConvertible = on", whereas thereactions with "isConvertible = off" are excluded from the mathematical model.
The parsedae module recognizes the variables prior to conversion of the reactions with "isConvertible = on". The
attribute of "totalAmount" of the element <specieAttribute> is introduced to determine whether the total
concentration for the species is constant (independent variable) or variable (time-dependent variable). When
"totalAmount" of protein is "variable", the protein is synthesized and decomposed according to TT, i.e., the total
concentration varies with time. When the total concentration of the protein with "totalAmount = constant" isindependent of time, protein synthesis and degradation are assumed not to occur. In addition, the parsedae module
distinguishes the species by the attribute "compartment" in the element of <specieAttribute> that indicates

which compartment the species is located.
At the first stage that generates ordinary differential equations, such as TT, CMA, GMA, and MM, the combined

attributes of "layer" in the elements of <listOfReactions>, "reactionType" and "groupID" of the elements of<Reactions>, and "modifierType" in the element of <listOfModifier> produce the appropriate differential
equations. The attribute of "layer", which has three values, "metabolic", "protein", and "gene", determines thebasic mathematical formula. The metabolic layer employs GMA or MM, the protein layer, which indicatessignal transduction pathways through protein interactions, uses CMA, and the gene layer uses TT. The
attribute of "reactionType", which consists of the various kinds of regulator-reaction equations, determines
which type of differential equations is employed in each type of reaction. The element of <listOfModifier>presents concrete regulator information, where the attribute of "modifierType" indicates the type of modifiers,"enzyme", "activator", or "inhibitor". The value of "enzyme" indicates general enzyme reactions, whereas"activator" or "inhibitor" is used for "reactionType = transcription". In transcription, an activator binds to an
enhancer to enhance the transcription of a gene, while an inhibitor binds to a promoter to repress the transcription ofa gene. The attribute of "groupID" is required to group a series of reactions, where the multiple regulators with thesame "groupID" act on the identical reaction. They are usually used for the "gene1" transcription. For example,

when multiple transcription regulation factors (TRFs) act on the transcription of a gene, theseregulator-reaction equations with the same "groupID" generate one transcription regulation equation withmultiple TRFs.
In order to express the spontaneous degradation of proteins and mRNAs in vivo, the attribute of"decomposition" is employed in the element of <specieAttribute>. Generally, proteins and RNAs
constitutively degrade in vivo due to factors including proteases ("decomposition =on "), whereas DNAs do
not degrade ("decomposition =off "). When the species consists of multiple components, the attribute of
12
KIT, Kurata Laboratory

"decomposition" has to be searched recursively, while checking whether "decomposition" of each component is on
or off.
In the second process, the attributes of "specieClass" "massBalance", and "totalAmount" in the element of<specieAttribute> play intrinsic roles in converting the ordinary differential equations into DAEs. The attributeof "specieClass" indicates the kinds of species, such as DNAs, RNAs, proteins, complexes, and metabolites.
This attribute of "specieClass" determines which components are needed to make the mass balance equations
for each species, e.g., "monomer", "modified", "enhancer", or "promoter". Binding multiple molecules without
any stoichiometric change generates "complex", whose components are listed in <listOfElements>. The
parsedae module automatically names the complex by joining its associated elements with a colon (:), e.g., PI:PII.
The value of "modified" indicates the molecules that are chemically modified through modification processes suchas phosphorylation and acetylation, which are accompanied with stoichiometric changes. The value of "modified" isautomatically named by joining the modifying molecule to the modified one using a hyphen (-), e.g.,
phosphorylated NRI is expressed as NRI-P. For example, the molecules of PI:PII, NRI-P:Enhancer, and
NRI-P:NRII are "complex", whereas NRI-P and NRII-P are "modified". The element of "modifierComplex" meansan active complex consisting of an enzyme and substrates, which is automatically generated when an enzyme
reaction is defined. The attribute of "massBalance" determines whether total mass balance equations are madefor the species with "massBalance = on", e.g., monomer, modified, enhancer, or promoter. TPP applies thequasi-steady state or rapid equilibrium approximation to the differential equations for the species with"specieClass = complex". For the species with "massBalance = on", the mass balance equations are generated
that sums all the complexes containing the species as components. The total concentration of the species isnamed by adding "T" to the head of species' name, whereas species' name without "T" indicates a free molecule.
In order to combine the gene and protein layers with the metabolic layer, the attribute of "isTotal" is presented todetermine if a modifier acts as a free species or as the total amount of the species. It is a critical attribute when theconcentration of the enzyme is defined as a time-dependent variable in the metabolic layer. With "isTotal = true" foran enzyme, the total concentration of the enzyme is employed for GMA or MM type. Conversely, with "isTotal =
false" for an enzyme, the free concentration of the enzyme is used for the metabolic layers. This rule enables one tocombine the three layers, i.e., metabolic and gene regulatory networks.

3 Conversion to mathematical models
3.1 TT (ordinary Transcription and Translation)
In the gene layer, transcription and translation occur in complicated manners, which involve a large number ofnucleic acids, amino acids, RNAs and proteins. It is not practical to describe all possible reactions, and it is difficult
to convert the transcription or translation (gene layer) into CMA based on their detailed molecular mechanism.
Thus, we apply a general type of differential equations suitable for transcription and translation equations. For
transcription, we conveniently divide two types: gene0 and gene1. For gene0, no transcription regulation factor
involves transcription; for gene1, activators or inhibitors regulate the transcription.

13
KIT, Kurata Laboratory3.1.1 Transcription without regulators (gene0)
Assuming that no gene regulation involves transcription, the transcription occurs at a constant rate, which isprovided by:Regulator-reaction equations: gene( A) → mRNA( A)(A1)d [mRNA( A)]= km [ gene( A)]
dt

(A2)

Differential equations:where km is the transcription rate constant.
3.1.2 Transcription with regulators (gene1)
For gene1, since activators or inhibitors regulate the transcription, regulator-reaction equations are provided by:

Activatori : Enhanceri− >> Gene → mRNA .

(A3)

Suppresori : Promoteri− Gene → mRNA(A4)

The synthesis rates are converted into:nmnm[ Suppressorj : Promoterj ][ Activatori : Enhanceri ]d [mRNA](+or×)∑ or ∏ (1 −)} ⋅ [Gene]= km ⋅ {∑ or ∏dt[TEnhanceri ][TPromoterj ]
i =1
j =1
i =1
j =1

(A5)

d [ Protein]= k p ⋅ [mRNA]
dt

(A6)

where km is the transcription rate constant, and k p is the translation rate constant. Complicated transcriptions
can be described by using the GUI editors.
3.1.3 Translation
For example, the protein A is synthesized from mRNA(A), which is given by:

d [ A]= k p [mRNA( A)] ,
dt

(A7)

where k p is the translation rate constant.

3.2 CMA (Conventional Mass Action)
Regulator-reaction equations can be divided into binding reactions and conversion reactions.
3.2.1 Binding reactions
For example, a binding reaction is provided by:
14
KIT, Kurata Laboratory→
α1 A1 +α2 A2 + K +αn An ←α1 A1 :α2 A2 : K :αn An ,

ka

(A8)

ks

where α i is the stoichiometry.
nd [ Ai ]α= α i − ka∏ [ A j ] j + kd [α1 A1 :α2 A2 : K :αn An ]
dt
j =1



(i = 1～n)nd [α1 A1 :α2 A2 : K :αn An ]α= ka∏ [ A j ] j − kd [α1 A1 :α2 A2 : K :αn An ] ,
dt
j =1

(A9)

(A10)

where ka is the binding association rate constant, and kd is the dissociation rate constant. The amount of thecorresponding species on the left and right sides has to be the same. The parsedae checks the stoichiometry ofreactions. In addition, the parsedae module checks the components of the complex. For the binding reaction
equations:→→
A1 + ( A2 : A3 ) ←A1 : ( A2 : A3 ) and A1 + ( A2 : A3 ) ←A1 : A2 : A3


ka

ka

ks

ks

(A11)

are allowed, but the equation:
→
A1 + A2 + A3 ←A1 : ( A2 : A3 )

ka

(A12)

ks

is not allowed, when the parsedae outputs the error.
3.2.2 Conversion reaction
The reaction:kx
α1 A1 : α 2 A2 : K : α n An →
β1B１ + β 2 B2 + K + β m Bm

(A13)

is converted into CMA:

d [ Bi ]= β i kx[α1 A1 :α2 A2 : K :αn An ] (i = 1～n)dt
d [α1 A1 :α2 A2 : K :αn An ]= − kx[α1 A1 :α2 A2 : K :αn An ] ,dtwhere kx is the conversion rate constant.

(A14)
(A15)

3.2.3 General enzyme reaction
General enzyme reactions except for transcription and translation can be provided by:ka
→ α A :α A : … : α 
α1 A1 +α2 A2 + … +αn An ←→ β1 B１ + β 2 B2 + … + β m Bm , (A16)
1 1
2 2
n
kx
kd

15
KIT, Kurata Laboratorywhere Ai is the substrates, Bi the products, α i andβ i are the stoichiometric coefficients, ka is theassociation rate constant and kd is the dissociation rate constant, and k x is the conversion rate constant. This

reaction is converted into CMA as follows:nd [ Ai ]α= α i  −ka ⋅ ∏ [ Aj ] j + kd [α1 A1 :α2 A2 : … : αn An ]
dt
j =1



(i = 1,..., n) .

(A17)

n
d [α1 A1 :α2 A2 : … : αn An ]α= ka∏ [ Aj ] j − kd [α1 A1 :α2 A2 : … : αβi An ] − βi kx[α1 A1 :α2 A2 : … : αn An ]
dt
j =1

(A18)

d [ Bi ]= β i ⋅ kx[α1 A1 :α2 A2 : … : αn An ]dt

(i = 1,..., m)(A19)

3.2.4 Transport
The species A is transported from one compartment to another compartment (A’)

d [ A]= − ktr[ A] ,
dt

(A20)

d [ A' ]= + ktr[ A]
dt

(A21)

where ktr is the transport rate constant. Active transport with an enzyme is processed in the same manner as aconversion with an enzyme.
3.2.5 Spontaneous decomposition
Proteins and RNAs are generally degraded in vivo. Thus, we omit the regulator-reaction equations for spontaneous
decomposition. Instead, we prepare the attribute of decomposition for species. With "decomposition = on" for a
species, the parsedae module adds the spontaneous degradation term to the differential equation for the species, as
follows:

d [ A]= − kpd [ A] ,
dt

(A22)

where kpd is the degradation rate constant. For mRNA, kmd is defined as the degradation rate constant for mRNAs.
Notice that the reaction type of "degradation" is different from the spontaneous decomposition.
When the components have complexes, the parsedae checks the attributes of decomposition recursively, while
checking whether "decomposition" is "on" or "off". For example, when the decomposition of A, B, C, A:B, B:2D,(A:B):(2(B:2D)):(3(C:D)) is on, that of D, C:D is off, the differential equations for the respective components
include the following terms :i)d[A]/dt= - kpd[1] * [A]ii)d[B]/dt= - kpd[2] * [B]16
KIT, Kurata Laboratoryiii) d[C]/dt= - kpd[3] * [C]iv) d[D]/dt= + 2*kpd[5] * [B:2D] + 4*kpd[6] * [(A:B):(2(B:2D)):(3(C:D))]v)

d[A:B]/dt= - kpd[4] * [A:B]vi) d[B:2D]/dt= - kpd[5] * [B:2D]vii) d[C:D]/dt= + 3 * kpd[6] * [(A:B):(2(B:2D)):(3(C:D))]viii) d[(A:B):(2(B:2D)):(3C:D)]/dt= - kpd[6] * [(A:B):(2(B:2D)):(3(C:D))]

Since the decomposition for the species of (A:B):(2(B:2D)):(3(C:D)) is on, the decomposition terms are added to
the respective differential equation (vii). Next, the decomposition of its components of A:B, B:2D, and C:D ischecked. Since the decomposition of A:B, A, and B is on, no source term is added to the equations (i,ii,v). The
decomposition of B:2D is on, the decomposition term is not added to the respective differential equation (iv).
However, since the decomposition of D is off, and the decomposition of B:2D releases 2*D, the source term for D
is added to the equation (iv). Since C:D is not decomposed despite the decomposition of (A:B):(2(B:2D)):(3(C:D)),the source term for C:D is added to the equation (vii).

3.3 GMA (General Mass Action)
Generally, metabolic reactions can be described using General Mass Action (GMA) or the Michaelis-Menten
equations (MM). GMA assigns each reaction to a mathematical term having a kinetic rate constant and power
coefficients regarding the concentrations of metabolites and enzymes, as described by:nd[ X i ] mf= ∑ kik ∏ [ X j ] ijk
dt
k =1
j =1

(A23)

where X1~Xn are the dependent variables, whose values vary with time, Xn+1~Xn+m are the independent variables,
whose values are fixed as constants, and kik is the rate constant. We show an example for mathematical conversion
for metabolic reactions provided by:

E −o S1 + S 2 > P1 + P2 .

(A24)

The GMA conversion cancels the concentration of the enzyme-metabolite complex of E:S1:S2, thus the differential
equation for E:S1:S2 is not produced, as follows.

dS1= − kxg1 E f1 S1 f 2 S2 f3
dt(A25),dS2= − kxg1 E f1 S1 f 2 S2 f3
dt(A27),

dP1= + kxg1 E f1 S1 f 2 S 2 f3
dt

(A28)

dP2
= + kxg1 E f1 S1 f 2 S2 f3
dt(A29),

where kxg1 is the kinetic rate constant, fi are the power coefficients. If necessary, GMA produces the term for
spontaneous decomposition as follows:
17
KIT, Kurata LaboratorydS= − kxgS f .
dt

(A30)

3.4 MM (simplified Michaelis-Menten equations)
Generally, the Michaelis-Menten equations are able to describe complicated reactions with activators or inhibitors
based on their mechanism, where the concentrations of the enzyme-metabolite complexes are cancelled. TheCADLIVE Simulator further simplifies the Michaelis-Menten equations. This simplification enables one to convert
complicated metabolisms into simple forms efficiently. For example, the metabolic equations can be converted into
the simplified Michaelis-Menten type differential equations. The regulator-reaction equation:

E −o S1 + S 2 > P1 + P2

(A31)

is converted to :

dS1
S1
S2
= − Q1 E
dt
Kmich1 + S1 Kmich2 + S 2(A32),dS2
S1
S2= − Q1 E
dt
Kmich1 + S1 Kmich2 + S2

(A33),dP1
S1
S2= + Q1 E
dt
Kmich1 + S1 Kmich2 + S2

(A34),dP2
S1
S2= + Q1 E
dt
Kmich1 + S1 Kmich2 + S2

(A35),

where K mich1 and K mich2 are the Michaelis-Menten constants. If necessary, this conversion also produces theterm for spontaneous decomposition as follows,

dS
S
=−Q
dt
Kmich + S

(A36).

3.5 DAEs
In order to solve the stiffness in CMA expressing a gene regulatory network, the two-phase partition method (TPP)
has been proposed that divide biochemical reactions into two phases, binding and reaction phases, by applying
rapid equilibrium approximation or quasi-steady state approximation to complex formation, which greatly reduces
the number of kinetic parameters and remarkably accelerates the calculation speed. TPP is a powerful method thatis indispensable for mathematical simulation of gene regulatory networks.
3.5.1 Steady-state approximation I
The two-phase partition method begins with CMA and TT.
(1) The differential equation for the complexes that have been generated through a binding reaction is set as zero,
resulting in an algebraic equation. The generated equations are simplified. For example, when a + b + c + d = 0 andb – d = 0 is generated, the former is replaced by a + c = 0.
(2) For the species that do not involve binding reactions and those whose attribute of "massBalance" is off, theirdifferential equations do not change.
18
KIT, Kurata Laboratory

(3) For the monomers that involve binding reactions and whose attribute of "massBalance" are on, theconcentrations of the monomer and its derivative complexes are all added, producing the differential equations fortheir total concentration. For the species whose attribute of "massBalance" is off, they are not. For example, for the
species regarding the monomer A, the equations are made as follows:
dy[TA]/dt =differential equations for monomer A+ Σ(differential equations for the complex that contains A as a component x the number of A thatthe complex contains)(4) The mass balance equations for the species whose "specieClass" is "monomer" is provided by:[TA] = [A] +Σ([complex that contains A as a component] x the number of A that the complex contains)(5) After simplifying the DAEs, the differential equations whose right hand sides become zero are omitted, and thespecies are registered as constant values.
The species are classified as follows:
A) Constant
The species whose attribute of "totalAmount" is constant, and whose "massBalance" is off.
The species whose differential equation is zero (dy/dt=0).
B) Variables of algebraic equations
The species whose class is "complex".
The species whose class is "monomer" and whose "massBalance" is on.
C) Variables of differential equations
The total concentration of the species with the attribute of monomer and modified. (excludes the species whose
differential equation is zero (dy/dt=0).)3.5.2 Steady-state approximation II
The algebraic equations at the steady state are further simplified by setting kpd = 0.
3.5.3 Rapid equilibrium approximation
The algebraic equations at the steady state are further simplified by setting kpd, kx << ka, kd. In addition, theassociation constant of Kb is substituted for ka/kd. For example, the regulator-reaction equation:
mA + nB<-> mA:nB(A37)

is converted to the algebraic equation:[mA:nB] = Kb^(m*n)*[A]^m*[B]^n.

(A38)

For example, the general enzyme reaction (A16) is converted into DAEs using the rapid equilibrium

approximation, as follows,
19
KIT, Kurata Laboratoryn0 = − K ⋅ ∏ [ Aj ] j + [α1 A1 :α2 A2 : … : αn An ](A39),[TAi ] = [Ai ] +αi [α1 A1 :α2 A2 : … : αn An ]

(A40)

α

j =1

d [TAi ]= − β i ⋅ kx ⋅ [α1 A1 :α2 A2 : … : αn An ]dtd [ Bi ]= β i ⋅ kx ⋅ [α1 A1 :α2 A2 : … : αn An ]dt

(i = 1,..., n)(i = 1,..., m)(i = 1,..., m)

(A41)
(A42)

where K is the association constant that is provided by ( ka / kd ).
3.5.4 Complexes in TPP
We illustrate in detail how TPP derives the differential equations for the total amount of specific species and theirmass balance equations. The regulator-reaction equation:
A + 2B <-> A:2B

(A43)

is converted into the CMA:
d[A]/dt = - ka[A][B]^2 + kd[A:2B](A44)

d[B]/dt = - 2ka[A][B]^2 + 2kd[A:2B](A45)

d[A:2B]/dt = + ka[A][B]^2 – kd[A:2B](A46).

We apply TPP to the CMA by setting the equation (A46) as zero, producing the differential equations for the total
concentrations (TA, TB) for the components of A and B. In order to make the differential equation for TB, werecursively search all the components of the complexes to find the species that contain B as a component, and addthem to the differential equation. They are multiplied by the number of B contained in the complexesd[TB]/dt = the right hand side of (A45) + the right hand side of (A46) x 2= Σ(the right hand sides of the CMA for the species that contains B as a component x the number of B
that the species contains)
In the same manner, the mass balance equation for TB is produced as follows:[TB] = [B] + [A:2B] x 2= Σ(the concentrations of the species that contain B as a component x the number of B that the species
contains)

3.6 S-system
Sensitivity and stability analysis shows how a biochemical system responds to perturbations or uncertainties.
S-system is able to analyze the various sensitivities of the system at the steady state in symbolic form. In thesimulator, ordinary differential equations such as TT, CMA, GMA, and MM can be converted into S-system toanalyze the sensitivities and stability at the steady state.
The ordinary differential equations are divided into the positive terms and negative terms:
20
KIT, Kurata Laboratory

dX i= Vi + − Vi − ,dt(i = 1,....,n,....n+m)（A47）,

where Vi+ is the sum of positive terms, and Vi− is the sum of negative terms. Generally, S-system is given by:
n+m
n+m
dX i
h
g
=αi ∏ X j ij −βi ∏ X j ij
dt
j =1
j =1

（A48）,

where X1 - Xn are the dependent variables, whose values vary with time, and Xn+1 - Xn+m are the independent
variables, whose values are fixed as constants. When the concentrations of the dependent variables are given at thesteady state, the coefficients ( gij , hij , α i , β i ) of S-system are solved from Eq. (A48) in symbolic form:

gij =

hij =

∂Vi + X j
∂X j Vi +∂Vi − X j
∂X j Vi −V+
αi = n + m ig∏ X j ij(A49),(A50),(A51),j =1

V−
βi = n + m ih∏ X j ij(A52),

j =1

where gij and hij are the kinetic order and α i andβ i are the rate constants.

3.7 Connection between layers
When DAEs express the gene-protein layers, and GMA or MM expresses the metabolic layer, dependent variables
for the regulators that appear as modifiers in the metabolic layer correspond to dependent variables for thegene-protein layers. Table A4 shows how to connect the corresponding variables between the gene-protein layer
and the metabolic layer. The combination of the values of the attributes of "isToal", "massBalance" and the method
for conversion determines the variables that appear as a modifier in the metabolic layer.

21
KIT, Kurata Laboratory

Table A4 Connection between layers
Conversion method for
massBalance isTotal The variable of E that appears as a modifier
Gene-Protein layer
in the metabolic layer
true
Total concentration of enzyme: T_E.cyt
onfalse
y[E.cyt]
CMA
true
Total concentration of enzyme T_E.cyt
off
false
y[E.cyt]true
y[TE.cyt]onfalse
x[E.cyt]
DAEs (TPP)
true
Total concentration of enzyme T_E.cyt
off
false
y[E.cyt]
The total concentration of the enzyme E is the total amount of E and the complexes that contain E as acomponent. For example, when TE consist of E, E:A:B, and E:C, the total concentration for the enzyme E is givenby:T_E.cyt = y[E.cyt] + y[E:A:B.cyt] + y[E:C.cyt](A53)

4. Selection of conversion methods
Selection of conversion methods is carried out for the two systems: gene-protein layer and metabolic layer.

4.1 Gene-protein layer
For the gene-protein layers, users select four types of conversion methods: CMA, TPP_STEADYSTATE_1,
TPP_STEADYSTATE_2, and TPP_RAPID, which are suitable for converting the regulator-reaction equations thatgenerate various complexes. The gene-protein layer can be converted into GMA or MM, but it is not practicallyuseful. Thus, such conversions are error.

4.2 Metabolic layer
For the metabolic layer, GMA and MM are generally employed. It is possible to convert the metabolic layer byusing CMA or TPP, but it is not useful practically.

4.3 Connection between the layers
Since the different conversion methods are applied to the gene-protein and metabolic layers, the differential
equations for the species common to both layers have to be combined to make one differential equation, whichenables one to simulate the three layers together.

22
KIT, Kurata Laboratory

B. CHECKDAE
1 Objectives
The variables employed in the checkdae file are named by the index made of the combination of species' name andits localization, as exemplified by y[RNAP.cyt], so that users can understand the meaning of the variables. However,in order for computers to simulate mathematical equations (TT, CMA, GMA, MM, DAEs, S-system), it is required
to replace the indexes of variable names by sequential integers, to make the user functions for mathematicalequations, and to produce their Jacobian function that is the partial derivatives of the user functions with respect toeach variable in symbolic form. The checkdae module generates computer-readable parameters, mathematical
functions, and their Jacobian. These functions and parameter files are described in the C programming language.
The checkdae module examines if a checkdae file (mathematical equations), which the parsedae module has
generated or users have edited manually, is acceptable for converting it into simulation programs (user functions
and parameters). The checkdae module generates not only the user functions for mathematical equations and theirassociated Jacobian, but also temporary mathematical expressions such as fluxes. In addition, the checkdae module
is able to convert ordinary differential equations (TT, CMA, GMA, MM) into S-system at the steady state,
generating the user functions for S-system with its associated coefficient functions and their parameters.
Simultaneously, it creates the parameter files for setting initial values and kinetic parameters.
In summary, from a checkdae file, the checkdae module creates the five files: a parameter file, a user function file,a temporary mathematical expression file, an S-system coefficient file, and an S-system user function file.

2. Function
Finishing checking the grammar of a mathematical model (checkdae file), the checkdae module divides it into thenodes consisting of one operator in order to expand them on a binary tree structure. The checkdae first examines if
the employed parameters and variables are defined in a checkdae file. Second, it converts the index of species'names to the numerical one, and solves the partial derivatives of the mathematical equations to obtain their Jacobianin symbolic form, resulting in a user function file with a parameter file that contains initial values and kinetic
constants. The use of GMA or MM conversion creates another file that contains the temporary mathematicalexpressions for its fluxes. When S-system is employed, it produces an S-system coefficient file and an S-systemuser function file.

2.1 Grammar check for mathematical models
The checkdae examines the correctness of the checkdae file (mathematical equations) and temporary mathematicalexpressions from the standpoint of grammar, which are restricted to five operators consisting of addition (+),
subtraction (-), multiplication (*), division (/), and power (^). The symbols except numbers are all regarded as
variables.

23
KIT, Kurata Laboratory2.2 Expansion to mathematical equation trees
Mathematical equations are expanded on binary trees whose nodes consist of an operator. The subsequent
operations are carried out based on the binary trees.

2.3 Function for checking parameters and variables
This function checks if the variables, parameters, and temporary mathematical expressions employed inmathematical equations have been defined in the checkdae file.

2.4 Function for expanding temporary expressions
The temporary mathematical equations are expanded to generate a user function, and then the temporary equations
are removed from the user function.

2.5 Function for indexing variables and parameters
The sequential integers are provided for all the variables and parameters that users have defined.

2.6 Function for replacing the indexes
The numerical indexes are substituted for the variables named with the combination of species' name and its
compartment.

2.7 Symbolic partial differentiation
The partial derivatives are solved with respect to the variables of interest in symbolic form.

2.8 Function for S-system conversion
Ordinary differential equations are converted into S-system. The detailed explanation is described elsewhere.

2.9 Parameter survey
Users are able to explore the parameter spaces by using the parameter file, where the values of parameters arevaried in the arithmetic or geometric series. Details are described in the parameter file (3.2.2).

3 Specification of input and output
3.1 Specification of input
The checkdae module reads the checkdae file (mathematical equations) to make the user functions and theirassociated parameter files necessary for simulation. The details of the checkdae file are explained in the instructionof the parsedae module.

3.2 Specification of output
3.2.1 User function file
The user function consists of mathematical equations and their Jacobian. The Jacobian is obtained in symbolic form
24
KIT, Kurata Laboratory

by solving the partial derivatives of mathematical equations with respect to each variable. The user function
contains two functions: usr_fvec and usr_fjac, as follows.
/** TITLE:test* INFO:test* CONVERSION TYPE*
GENE-PROTEIN:TPP_RAPID*METABOLIC:NONE*/#include "life.h"void usr_fvec(double fvec[], double y[], double Gene, Parameter *p){
double *constantPlayer= p->val[1];double *Kb= p->val[2];double *kp= p->val[3];double *kpd
= p->val[4];double *km= p->val[5];double *kmd
= p->val[6];
fvec[
fvec[
fvec[

1] = y[12] - (y[1] + y[6] + y[7] + pow(y[8], 1.5));2] = constantPlayer[2] - (y[2] + pow(y[6], 2) + pow(y[7], 2) + pow(y[8], 3));3] = constantPlayer[3] - (y[3] + pow(y[7], 1.5) + y[10]);--- omission --fvec[ 12] = -kpd[1]*y[1] - kpd[2]*y[6] - kpd[3]*y[7] - 1.5*kpd[4]*y[8];fvec[ 13] = kp[1]*y[11] - kpd[5]*y[5] - kpd[6]*y[9] - kpd[7]*y[10];}void usr_fjac(double fjac[][MAX_VAR_NUM], double y[], double Gene, Parameter *p){
double *constantPlayer= p->val[1];double *Kb= p->val[2];double *kp= p->val[3];double *kpd
= p->val[4];double *km= p->val[5];double *kmd
= p->val[6];
fjac[
fjac[

1][
1][

1] = -1;6] = -1;--- omission --fjac[ 13][ 9] = -kpd[6];fjac[ 13][ 10] = -kpd[7];fjac[ 13][ 11] = kp[1];}Equations are indexed in the sequential order of algebraic and differential equations. Variables are marked by "y".
The fjac is described only when its value is nonzero.
3.2.2 Parameter file
The checkdae module outputs the parameter file for setting initial values and kinetic parameters. The line up of
PARAM_NAME corresponds to the numerical index order of the parameters in the user function. When users gave
25
KIT, Kurata Laboratoryvalues to parameters and initial concentrations in the chekdae file, the values appeared at the corresponding order of
PARAM_NAME.
##### Model ################# TITLE:test# INFO:test# CONVERSION TYPE#
GENE-PROTEIN:TPP_RAPID#METABOLIC:NONE
N_VAR; 13; # num of variables(all)
N_ALGEBR; 10; # num of variables(Algebraic Eq)#Y_START ; index; initial_value; tag #comment
Y_START ;1; 1.2000e+00; A.cyt#
Y_START ;2; 2.1000e+00; B.cyt#
Y_START ;3; 3.1000e+00; C(pro).cyt--- omission --Y_START ; 11; 1.0000e-01; mRNA(test).cyt
Y_START ; 12; 1.7650e+00; TA.cyt #Y_START ; 13; 6.0000e-02; Ttest.cyt ####PARAM;name;index;val;val_start;num_survey;D/R/S;change-val;GA_start;GA_end;tag #comment
PARAM;constantPlayer;
1; 1.2000e+00;;0;D;;;; G(test).cyt#
PARAM;constantPlayer;
2; 3.2300e+00;;0;D;;;; TB.cyt #--- omission --PARAM
;kmd
;
1;
;;0;D;;;; decomposition_rate_constant_mRNA(test).cyt##T_EVENT ; name; index; time; value; #comment######## Don't edit the following lines. ######PARAM_NAME ;constantPlayer #PARAM_NAME ;Kb#
PARAM_NAME ;kp#
PARAM_NAME ;kpd#
PARAM_NAME ;km#
PARAM_NAME ;kmd#

The checkdae module also implements the function for surveying kinetic parameters to investigate the behaviorof a model over a specific range of the values of some parameters. This allows users to effectively study thedependency of the model's behavior on the kinetic parameters, and is therefore a very effective means of forecasting
the effects of parameter perturbations.
Details of the #PARAM line are explained as follows,
name: the name of the parameter
index: the index of the parameter
val: the value of the parameter
val_start: the starting value of the parameters, which must be set as the same value as "val" when the parameter
survey is carried out.
26
KIT, Kurata Laboratorynum_survey: the number of surveyed parameters. When the parameter survey is not performed, set one (1).
D/R/S: The option for the parameter survey. "D" indicates arithmetic series, "R" indicates geometric series, and "S"means the parameter search by GAs or RS (Details are written elsewhere).
change_val: changing value.
GA_start; The starting value of the explored variable region, which is employed for search for GAs or RS.
GA_end; The ending value of the explored variable region, which is employed for search for GAs or RS.
For example, when users change the parameter kx[3] as 1, 2, 3, and 4 in the arithmetic series, the line of #PARAM
is provided by:
PARAM

;kx;

3;

1;

1;

4;

D;

1;

;

;

tag #

When users change the parameter kx[3] as 1, 2, 4, and 8 in the geometric series, the line of #PARAM is provided
by:
PARAM

;kx;

3;

1;

1;

4;

R;

2;

;

;

tag #

3.2.3 Temporary mathematical expression file
When temporary mathematical equations such as a flux and total enzyme are defined in the checkdae file, thecorresponding user function is generated automatically. The temporary expressions appear in the checkdae file,when GMA or MM conversion is selected. Mathematical equations for fluxes are shown as follows.
/** TITLE:test* INFO:test* CONVERSION TYPE*
GENE-PROTEIN:TPP_RAPID*METABOLIC:NONE*/#include "life.h"void usr_flux(double total[], double flux[], double y[], double Gene, Parameter *p){
double *constantPlayer= p->val[1];double *Kb= p->val[2];double *kp= p->val[3];double *kpd
= p->val[4];double *km= p->val[5];double *kmd
= p->val[6];double *Kg= p->val[7];double *f = p->val[8];flux[1] = Kg[1]*pow(y[1], f[1]); /* flux_A.cyt_to_B.cyt *//* No Total expression. */}3.2.4 S-system coefficient file
The S-system coefficient file is used for calculating the coefficients of S-system at the steady state, which are sent
27
KIT, Kurata Laboratory

to the program for the sensitivity and stability analysis of S-system.
/** TITLE:Glycolysis.plc* INFO:Glycolytic-Glycogenolytic Pathway Model in Chapter 11 PLUS*/
#include "life.h"void get_GHdata(double **gd, double **gi, double **hd, double **hi,double *v_plus, double *v_minus, double *y, Parameter *p){
double *constantPlayer= p->val[1];v_plus [ 1] = 0.0778843*pow(constantPlayer[1], 0.66)*constantPlayer[3];gi[ 1][ 1] = 0.66*pow(constantPlayer[1], -0.34)*0.0778843*constantPlayer[3];gi[ 1][ 3] = 0.0778843*pow(constantPlayer[1], 0.66);
v_minus[ 1] = 1.06271*pow(y[1], 1.53)*pow(y[2], (-0.59))*constantPlayer[4];hd[ 1][ 1] = 1.53*pow(y[1], 0.53)*1.06271*pow(y[2], (-0.59))*constantPlayer[4];hd[ 1][ 2] = (-0.59)*pow(y[2], (-0.59) - 1)*1.06271*pow(y[1], 1.53)*constantPlayer[4];hi[ 1][ 4] = 1.06271*pow(y[1], 1.53)*pow(y[2], (-0.59));--- omission --gi[ 3][ 5] = 0.000793456*pow(y[2], 3.97)*pow(y[3], (-3.06));v_minus[ 3] = 1.05881*pow(y[3], 0.3)*constantPlayer[6];hd[ 3][ 3] = 0.3*pow(y[3], -0.7)*1.05881*constantPlayer[6];hi[ 3][ 6] = 1.05881*pow(y[3], 0.3);}
void get_GHnum(int *n_var, int *n_const){*n_var = 3;*n_const = 7;}

4 Data structure
4.1 Main structures
EqSystem_T|+-int n_var
total number of variables = n_varX + n_varY+-int n_varX
number of x variables+-int n_varY
number of y variables+-int n_mid
number of temporary mathematical expressions+-int n_equation
number of equations (n_varX + n_varY)+-int n_param number of parameter species (containing constantPlayer)|+-Variable_T *var(n_var) variable arrays
|
|
|
+-int
index
28
KIT, Kurata Laboratory|+-char*tag
molecular name|+-double
val
|
+-int
f_isX
x or y flags|+-int
f_isUsedflags for use or non-use|+-Expression_T *mid(n_mid)array for temporary mathematical expressions
|
|
|
+-int
f_isXorT(X/Y or) Total/Flux flag|+-char*nameleft hand side equations (identification name)|+-char*expr
expanded equations|||+-Node_T *root
expanded equations in binary trees (list)||+-Expression_T *equation(n_varX + n_varY) array for mathematical expressions
|
|
|
+-int
f_isXorT
x or y flags|+-char*name
Y: it is used, X: it is not used.
|+-char*expr
expanded equations|||+-Node_T *root
expanded equations in binary trees (list)|+-ParamBase_T *param_base(n_param)
array for parameter species|+-char*labelka, kb, etc+-int
max_index
maximum index+-int
f_isUsedflags for use or non-use|+-Parameter_T *param(max_index)
array for each parameter species|+-char*tag
content+-doubleval+-int
f_isUsedflags for use or non-use

4.2 Node_T structure
In pdiff.h, the definition of the Node_T structure, which is key to mathematical derivations, is provided as follows.
typedef struct _node {
char
op;
int
rank;
int
f_bracket;
int
f_chunk;struct _node*parent;
struct _node*childL;char*charL;int
rank_charL;
struct _node*childR;char*charR;int
rank_charR;char*expr;
int
rank_expr;char*diff;
int
rank_diff;} Node_T;
29
KIT, Kurata Laboratory

This structure is generated with respect to an operator (+ , -, *, / , or ^). Mathematical equations are expressed bylinking the structures in a binary tree. Binominal operations are carried out, but it is possible to define '-' as a single
operator. In this case, neither childL nor charL is defined. In the other cases, operands are defined in both the rightand the left. For example, the left operand is defined either as charL or childL. If the operand is provided as charL,
rank_charL defines the value that identifies either numerical values (R_5VAL) or strings (R_4CHAR), wherechildL is NULL. Conversely, if the operator is defined as childL, charL is NULL and rank_charL is provided by
R_NULL. The right operand is defined in the same manner.
The others are as follows.
char

op;
define an operator of the node (+, -, *, / , or ^)

int

rank;
define the priority value of operators (R_1PLUS、R_MULT、R_3POW) to the above operator.

int

f_bracket;
frag necessary for handling "(" on the way to expansion of equations.

intf_chunk;flag showing that the operation of the node employs brackets.

struct _node*parent;
pointer to the parent
char*expr;
registering mathematical expressions

int
rank_expr;
showing the rank of the above mathematical expressions, which defines R_**.
char*diff;
registering mathematical expressions that are partially differentiated.

int

rank_diff;showing the rank of the above mathematical expressions, which defines R_**.

In order to identify the priority of an operator and the order of operations, rank** is provided as follows.
R_NULL
R_1PLUSR_2MULT
R_3POW
R_4CHAR
R_5VAL:not defined:addition or subtraction, or mathematical expression containing them
:multiplication or division, or mathematical expression containing them but excluding additionand subtraction:power, or mathematical expression consists of powers.
:variables (in mathematical expression):values (in mathematical expression)

5 Procedure
The relationships of functions are described in a tree structure. The functions followed by ‘…’ indicates that they
have been already called. The functions preceded by ‘…’ are recursive functions.
30
KIT, Kurata Laboratorymain+-my_init
setting file names+-Read_data
read data|+-split_dae_file
division of a checkdae file, read the numbers of data|+-allocate_area
memory allocation for data|+-read_param_base
read parameter labels.
||+-divide_line
Strings are sectioned by '=' or ';'|+-read_paramread parameters|| +-read_param_index
obtain the maximum number of parameters|||+-divide_line ...
|||+-Divide_variabledivide variable and parameter strings by '[' and ']'|||+-FindParamBaseAddr
examine if parameter labels are registered||+-read_param_data
read parameter data||+-divide_line ...
||+-Divide_variable ...
||+-FindParamBaseAddr ...
|+-read_var
read variables||+-divide_line ...
||+-Divide_variable ...
|+-read_mid_expression read temporary mathematical expressions|+-fix_equation_line
rearrangement of equations||+-Trim_tailtrim line ends||+-divide_line ...
|+-read_equation
read differential and algebraic equations|+-fix_equation_line ...
|+-divide_line ...
|+-Divide_variable ...
+-set_expression expand mathematical equations in a binary tree.
|+-Make_node
expand mathematical equations in a tree structure||+-create_node
create nodes||+-get_token
read mathematical equations||+-check_token
check the grammar of mathematical equations||+-add_emptyChildcreate blank child nodes|||+-create_node ...
||+-seek_bracket
search a node whose flag is '('|||+-... seek_bracket ...
||+-del_node
delete nodes||+-set_oprator
set operator information||+-insert_childNodeinsert child nodes|||+-create_node ...
||+-seek_same
search equivalent or less prior nodes|||+-... seek_same ...
||+-insert_parentNodeinsert a parent node|||+-create_node ...
||+-set_operanda set of operands|+-Set_expr
reconstruct mathematical equations based on a tree structure||+-... Set_expr ...
||+-del_bracket
delete brackets||+-add_bracket
add brackets|+-check_node
check nodes recursively|+-... check_node ...
|+-replace_variablereplace with numerical indexes||+-Divide_variable ...
||+-getVarIndex
obtain numerical indexes for variable names|+-check_detailcheck nodes|+-Divide_variable ...

31
KIT, Kurata Laboratory|+-findMidAddr
examine if temporary mathematical expressions are registered|+-expand_mid
expand temporary mathematical expressions||+-Make_node ...
||+-Seek_root
obtain the address of a root node|||+-... Seek_root ...
||+-Set_expr ...
|+-findVarAddr
examine if variables are registered|+-FindParamBaseAddr ...
|+-setUsedParam
Add flags to employed parameters+-set_param_flag examine if parameter labels are used or not+-set_variable_index
assign numerical indexes to variables|+-getDiffEqAddr
obtain the address for registering differential equations+-set_numerical_index
replace with numerical indexes|+-check_node ...
|+-Set_expr ...
+-Write_usrfunc
create user functions for simulator|+-write_header
write header information|+-write_param_label
write information that relates parameters|+-Make_diff
symbolic partial differentiation|+-... Make_diff ...
|+-diff_add
partial differentiation for addition||+-del_bracket ...
|+-diff_subtract
partial differentiation for subtraction||+-del_bracket ...
|+-diff_multiply
partial differentiation for multiplication||+-multiply
multiplication||+-del_bracket ...
|+-diff_divide
partial differentiation for division||+-multiply ...
|+-diff_power
partial differentiation for power||+-multiply ...
||+-del_bracket ...
|+-add_bracket ...
+-Write_paramFile
create parameter files|+-write_header ...
+-Write_flux
create temporary mathematical expression files|+-write_header ...
|+-write_param_label ...
|+-Set_Ssystemwrite S-system user function files+-write_GEhead
write the header part of the file that defines g and h|+-Write_header ...
|+-Write_param_label ...
+-make_Sdata
write the core part of S-system user functions|+-Set_expr ...
|+-separate_terms
divide differential equation into V+ and V|+-get_preGH
writhe g, h, V+, and V||+-Del_tree
delete unnecessary expansion trees|||+-... Del_tree ...
||+-Make_node ...
||+-Set_expr ...
||+-Make_diff ...
|+-write_Sfvec
write fvec in a temporary file|+-write_Sfjac
write fjac in a temporary file+-write_GEtail
write the end part of the file that defines g and h+-write_S_usrfunc
write S-system user function file+-Write_header ...

32
KIT, Kurata Laboratory

+-write_S_funcs

write S-system differential equations

Main procedures are explained as follows.

5.1 Input procedure (Read_data)
5.1.1 Reading checkdae files: data split and data number acquisition (split_dae_file)
A checkdae file is read as input data, and divided according to the identification lines. The divided data are sent to atemporary file. The numbers of various data are obtained by counting semicolons. Each valid data has one
semicolon.
5.1.2 Index check and memory allocation (allocate_area).
After checking the numbers the have been obtained by split_dae_file, the memory necessary for structure arrays forthe employed parameters are allocated.
5.1.3 Reading parameter labels (read_param_base)
After reading the parameter labels from the temporary file generated by split_dae_file, the parameter labels andconstantPlayer are registered in ParamBase_T.
5.1.4 Reading parameters (read_param)
The maximum number of each parameter species is read from the temporary file generated by split_dae_file, andthe memory necessary for the parameter arrays is allocated (by read_param_index). Then the data is read again andregistered in Parameter_T (by read_param_data).
5.1.5 Reading variables (read_var)
The variables are read from the temporary file generated by split_dae_file, and registered in Variable_T. The
correctness for the assigned indexes is checked.
5.1.6 Reading temporary mathematical expressions (read_mid)
The temporary mathematical expressions are read from the temporary file generated by split_dae_file, and remade
by fix_equation_line into another temporary file. The remade data are read, and registered in Expression_T.
5.1.7 Reading differential equations and algebraic equations (read_equation)
The files for differential equations and algebraic equations that have been generated by split_dae_file are rewritten
by fix_equation_line in other temporary files. From which, equation data are read and registered in Expression_T.
In processing algebraic equations, the right-hand equations are moved to the left. The resultant left-had side isregistered as expr. In processing differential equations, the right-hand side is registered in expr. The index of theleft-hand side equation is registered as name.

33
KIT, Kurata Laboratory5.2 Expansion of mathematical equations in a tree structure (set_expression、Make_node)
The temporary mathematical expressions, differential equations, and algebraic equations, which are registered as
expr in Expression_T, are expanded in a tree structure by Make_node, where Node_T is employed as nodes. In
Make_node, the root node is made as a current node. The strings of mathematical equations are sectioned by anoperator, an operand, or brackets (by the get_token), and their grammar is checked (by check_token). According to
the data, a new node is linked to the current node as a parent or a child. This process is iterated until the whole
equations are read. Consequently, expr registers the expanded mathematical equations in a node tree. After
expansion, mathematical expressions are reconstructed based on the generated tree (by Set_expr function).
Important functions in Make_node are follows.
5.2.1 Reading mathematical expression (get_token)
Strings of mathematical expressions are read by a character, grouped by an operator, an operand, or brackets, andregistered in the structure of Token_T. Token_T is as follows:
typedef struct _token {char
op;operator or parentheses: '+ ', '(', variable: ’a’, value: ‘v’
int
rank;
operator:R_PLUS、R_MULT or R_POW. others: R_NULL。
char*expr; pointers to strings of variables and values
int
bra;
counter for '(' (check_token)
int
cket;
counter for ')' (check_token)} Token_T;
The index within [] is recognized as characters. All the blanks are eliminated during pursing.
5.2.2 Checking mathematical expression (check_token)
The correctness of the data read by get_token and the balance of brackets are checked.
5.2.3 Adding a blank child node (add_emptyChild)
If the operator of the current node is set, a blank child node is added to the right. Otherwise, it is added to the left.
5.2.4 Inserting a child node (insert_childNode)A new child node is inserted as the right hand of the current node. The new node receives the right-handinformation of the parent node as the left-hand information.
5.2.5 Inserting a parent node (insert_parentNode)A new node is inserted as the parent node of the current node. The current node is the right-hand child node of the
new node.
5.2.6 Deleting nodes (del_node)
It the extra number of brackets are found, the unnecessary nodes are eliminated.
34
KIT, Kurata Laboratory5.2.7 Searching equivalent or less prior nodes (seek_same)
When the operator that has been read shows a high priority to the current node, it is defined as a child node of thecurrent node. However, when the operator shows a less priority or equivalence to the current, the nodes, whose
priority is equivalent or less than the operator, whose flag is '(', or whose operator is not defined, are searched in the
direction of the root, placing the operator node in the correct position. If the node is not found, return the root.
5.2.8 Reconstructing mathematical expressions (Set_expr)
Mathematical expressions are reconstructed recursively from the expanded binary tree by Make_node. Numericalvalues are calculated when they appear.

5.3 Checking variables (check_node、check_detail)All the variables, parameters, and temporary mathematical expressions that are defined in cahrL and charR of eachnode of a binary tree are recursively checked (by check_node). If there are non-defined values, the program willstop and display "error".

In the case of variables and parameters, the flag of f_isUsed is put up. In the case of

temporary mathematical expressions, they are expanded in a tree by expand_mid.

5.4 Expanding temporary mathematical expressions (expand_mid)
When a temporary mathematical expression is contained in a mathematical equation, the expression is expanded ina binary tree by Make_node. The generated tree is connected to the child node of the mathematical equation. Finally,
the bound expressions are reconstructed by Set_expr.

5.5 Assigning an index to a variable (set_variable_index)
Integer consecutive numbers are assigned to the employed variables.

5.6 Converting variables (set_numerical_index, check_node, replace_variable)
After searching the nodes that set variables (by check_node), the name-indexed variables are replaced by integerindexes (by replace_variable). After replacement, the resultant expressions are reconstructed by Set_expr.

5.7 Symbolic partial differentiation (Make_diff)
Each term for the expanded mathematical equations is partially differentiated from the root node. The process iscarried out recursively. The partial differentiation is applicable to the five operators:
addition:diff_add

subtraction:diff_subtract

multiplication:diff_multiply

division:diff_divide

power:diff_power
35
KIT, Kurata Laboratory

The Diff_T structure is employed to send and receive data.
typedef struct _diff {int
rank_exprL;
int
rank_exprR;
int
rank_diffL;
int
rank_diffR;
char
exprL[BUFFSIZE];
char
exprR[BUFFSIZE];
char
diffL[BUFFSIZE];
char
diffR[BUFFSIZE];} Diff_T;
Make_diff is called to make Jacobian, which is defined as a user function.

5.8 Output process
5.8.1 Output of user functions for simulators (Write_userfunc)
There are two types of user functions: usr_fvec and user_fjac. The usr_fvec consists of algebraic and differential
equations that have been expanded and replaced with numerical indexes. The user_fjac is the Jacobian, where thevariable is partially differentiated with respect to each variable. Note that the fjac whose value is zero is excluded.
First, the header part that has been obtained by split_dae_file is described as comments in the C language format,
and then usr_fvec and usr_fjac are output. In the user functions, the employed parameters are linked to theparameter structure (write_param_label). In user_fvec, expr in the Expression_T structure that stores mathematicalexpressions is output sequentially. In usr_fjac, Make_diff is applied to a differential equation with respect to eachvariable, and the resultant Jacobian is output in order.
5.8.2 Output of parameter files (Write_paramFile)First, the header part that has been obtained by split_dae_file is described as comments. Next, the initial value
setting part, parameter-setting part, and the data that relate them to the parameter labels are described. When theinitial values and parameter values are not set, blanks are output.
5.8.3 Output of temporary mathematical expression obtaining file

(Write_flux)

Instead of usr_fvec, usr_flux is output that is the function to defines temporary mathematical expression The total
concentration of enzymes is defined as the array "total", that of fluxes as array of "flux". The original temporary
mathematical expressions are output as comments on the right-hand side of the equations.

5.9 S-system conversion（Set_Ssystem）
The S-system conversion is carried out with respect to each differential equation by make_Sdata, as follows:

36
KIT, Kurata Laboratory

1. Differntial equations are redescribed by Set_expr.
2. All the terms of the right-hand side are divided into V+ and V-, respectively, by separate_terms.
3. V+ and V- are output by get_preGH. Made_diff generates and writes

∂V+
∂V−
and.
The partial∂X∂X

differential formula are described separately with respect to the dependent or independent variables.
4. write_Sfvec writes the differential equations of S-system with respsect to the dependent variables.
5. write_Sfjac generate and write the Jacobian.
6. The resultant expression are output as the S-system coefficient file and S-system user function file.

3.2.5 S-system user function file
The checkdae module converts ordinary differential equations into S-system, generating the user function of
S-system, as follows.
/** TITLE:Glycolysis.plc* INFO:Glycolytic-Glycogenolytic Pathway Model in Chapter 11 PLUS*//*
#include "life.h"void usr_fvec(double fvec[], double y[], double Gene, Parameter *p){
double *constantPlayer= p->val[1];double *kplus= p->val[2];double *kminus= p->val[3];double *g = p->val[4];double *h = p->val[5];
fvec[
fvec[

1] = kplus[1]*pow(constantPlayer[1],g[1])*pow(constantPlayer[3],g[2])- kminus[1]*pow(y[1],h[1])*pow(y[2],h[2])*pow(constantPlayer[4],h[3]);--- omission --3] = kplus[3]*pow(y[2],g[8])*pow(y[3],g[9])*pow(constantPlayer[5],g[10])- kminus[3]*pow(y[3],h[7])*pow(constantPlayer[6],h[8]);}void usr_fjac(double fjac[][MAX_VAR_NUM], double y[], double Gene, Parameter *p){
double *constantPlayer= p->val[1];double *kplus= p->val[2];double *kminus= p->val[3];double *g = p->val[4];double *h = p->val[5];
fjac[
fjac[
fjac[
fjac[

1][ 1] = - kminus[1]*h[1]*pow(y[1],h[1]-1)*pow(y[2],h[2])*pow(constantPlayer[4],h[3]);1][ 2] = - kminus[1]*pow(y[1],h[1])*h[2]*pow(y[2],h[2]-1)*pow(constantPlayer[4],h[3]);--- omission --2][ 3] = - kminus[2]*pow(y[2],h[4])*h[5]*pow(y[3],h[5]-1)*pow(constantPlayer[5],h[6]);3][ 2] = kplus[3]*g[8]*pow(y[2],g[8]-1)*pow(y[3],g[9])*pow(constantPlayer[5],g[10]);37
KIT, Kurata Laboratory

fjac[

3][3] = kplus[3]*pow(y[2],g[8])*g[9]*pow(y[3],g[9]-1)*pow(constantPlayer[5],g[10])- kminus[3]*h[7]*pow(y[3],h[7]-1)*pow(constantPlayer[6],h[8]);}38
KIT, Kurata Laboratory

C. SOLVER
1. Introduction
Two distinct types of simulation can be carried out: (i) time course simulation, where the values of variables are
determined as a time series; (ii) steady state analysis, where the values of variables are determined for a state inwhich metabolite concentration does not change. For differential equations, the CADLIVE Simulator implements
the Runge-Kutta method, the step-adaptive Runge-Kutta method, and the NDF that can be applied to highly stiff
differential equations. The NDF has been programmed based on ode15s of MATLAB. The Runge-Kutta or the
step-adaptive Runge-Kutta method is employed for non-stiff differential equations. The Newton-Raphson algorithm
is employed to solve algebraic equations. The combination of the Newton-Raphson algorithm with the Runge-Kutta
or with the step-adaptive Runge-Kutta method, or the NDF solves DAEs. The solver not only simulates dynamic
behaviors of biochemical networks, but also solves the steady state solutions by applying the quasi-steady state
approximation to all the differential equations, resulting in generating algebraic equations. The Newton-Raphson
method solves such algebraic equations to obtain the steady state concentrations. In order to enhance the calculation
rate, the CADLIVE Simulator employs the Message-Passing Interface (MPI) for parallel calculation.

2. Analysis types
As "Analsys type", users can choose either "Dynamic” or "Steady-state". "Dynamic" simulates the time evolutionof the concentrations by calculating DAEs, and "Steady-state" calculates the concentrations at steady state bysolving algebraic equations. The checkbox of "Parameter survey" determines if the simulator surveys the parameter
space. Checking the checkbox of "Parallel calculation" carries out parallel calculation that employs the Message
Passing Interface (MPI). Notice that the checkbox of "Parallel calculation" cannot be selected prior to checking theparameter survey.

3. Input for control data
The control parameters for the solvers are described as Table C1-C3. How to provide the value is shown in theexamples (CADLIVE web page).

39
KIT, Kurata Laboratory

Table C1 Control parameters for differential equations
Solver Type
Runge-Kutta
Runge-Kutta (Adaptive step-size)
NDF
Set time span and time step-size.
Start time
End time(Initial)time step-size
Monitoring interval

Constructed based on "Numerical Recipes inC" (Cambridge University Press)
Constructed based on ode15s of MATLAB
The start time of simulation
The end time of simulation
Initial time step size
The simulated time course is picked up by themonitoring interval value.

Table C2 Control parameters for algebraic equations
Solver type
Newton-Raphson Method
Maximum trial times
Tolerance for convergence of functions
Tolerance for convergence of variables
Ratio of changing parameterschange width calc. sensitivity (STD).

Constructed based on "Numerical Recipes inC" (Cambridge University Press)
The number of iterated calculations
Convergence condition
Convergence condition
When the value for kinetic parameter ischanged, the value should be varied by theratio of changing parameters. Because the greatchange in the kinetic parameter sometimes
causes the calculation to fail.
The ratio of perturbation that is given to
parameters is set as STD.
This is used for analyzing the sensitivity ofmodels.

Table C3 Other parameters
Other
G-value
Y default value

Always 1.0
default values for y (molecular concentrations)

40
KIT, Kurata Laboratory

D. MERGEPARAM
1. Objectives
It is laborious and time-consuming to set many parameters manually. The mergeParams module helps users setting
parameters by copying or merging the existing data, which greatly reduces such laborious parameter setting. The
mergeParams module updates the current parameter file by a new parameter file.

2. Function
The mergeParams module updates the values for the variables, constants, and parameters of the current parameter
file by a new one according to their tags (#comment). The tags are added to all the variables, constants andparameters in order to distinguish them. Comparing the new parameter values with the current parameter values, thecurrent parameters with the same tag as the new ones are updated, if necessary. When some values of the current
parameters have been already set, users can chose two types:
1) Only empty parameters are updated.
2) The values for all the current parameters are updated.

3. Input/output specification
The mergeParams module reads a new parameter file and the current file, and updates the current parameters by the
new ones. Both parameter files are text files, the lines to be merged are:・ Y_START: the line starting from "Y_START". the initial values of variables.

・ PARAM: the line starting from "PARAM".

the values of kinetic parameters.

This module processes only the lines to be merged, the other lines never be changed.

3.1 Y_STARTY_START ; index; initial_value; tag #comment
Each field is distinguished with semicolon";". The mergeParams skips the region preceded by "#". The tag"#comment" distinguishes the parameters to be merged, and "initial_value" is merged.

3.2 PARAM
PARAM;name;index;val;val_start;num_survey;D/R;change-val;GA_start;GA_end;tag #comment
Each field is distinguished with semicolon";". The mergeParams skips the region preceded by "#". The tag
41
KIT, Kurata Laboratory

"#comment" distinguishes the parameters to be merged, and "val" is merged.

4. Command
Users use the mergeParams program on command lines as follows.
mergeParams –IMathParam.txt –Rinput.txt [–B]
MathParam.txt is the current parameter file. The space between the file and the option is not allowed.
・ MathParam.txt

input/output : the current parameter file・ input.txt

input: the parameter file that updates the current one.

The flag –B indicates that only empty parameters are updated.
No flag indicates that all the current parameters are updated.

5. Error message"Too many parameters.""Invalid line:
LINE_TEXT""Cannot open FILENAME for reading.""Cannot read parameters. (file = FILENAME)""Cannot open FILENAME for writing.""Invalid argument !!""Too few arguments !!""Cannot import files !!""Cannot merge parameters !!""Cannot export parameters !!"42
KIT, Kurata Laboratory

E. SENSITIVITY AND STABILITY ANALYSIS BY S-SYSTEM
1. Objectives
Sensitivity analysis shows how a biochemical system responds to perturbations or uncertainties. Among manymethods for system analysis, the sensitivity and stability analysis is useful for characterizing the robustness ofmathematical models at a steady-state level, and allows one to determine which parameters have the most effect on
the dynamics, or which factors cause to the system to oscillate. S-system is employed to analyze the various
sensitivities and stability of the system at the steady state, because the use of S-system solves them in symbolicform. In this module, ordinary differential equations such as TT, CMA, GMA, and MM are converted into S-system
to analyze the sensitivities and stability.

2. Introduction of S-system
We show how S-system is derived from ordinary differential equations, and how the stability and sensitivity aredefined at the steady state.

2.1 S-system conversion
Ordinary differential equations are divided into the positive terms and negative terms:

dX i= Vi + − Vi − ,dt(i = 1,....,n,....n+m)(E1),

where Vi+ is the sum of positive terms, and Vi− is the sum of negative terms. Generally, S-system is given by:
n+m
n+m
dX i
g
h
=αi ∏ X j ij −βi ∏ X j ij
dt
j =1
j =1

(E2),

where X1,…, Xn are the dependent variables, whose values vary with time, and Xn+1,…,Xn+m are the independent
variables, whose values are fixed as constants. When the concentrations of the dependent variables are given at thesteady state, the coefficients (gij, hij, αi, βi) of S-system are solved from Eq. (E2) in symbolic form:

gij =

hij =

∂Vi + X j
∂X j Vi +∂Vi − X j
∂X j Vi −V+
αi = n + m ig∏ X j ij

(E3)

(E4)

(E5)

j =1

43
KIT, Kurata Laboratory

V−
βi = n + m ih∏ X j ij

(E6),

j =1

where gij and hij are the kinetic order and αi and βi are the rate constants. At the steady state:

dX i
= 0,
dt
we can take logarithms on both sides of Eq. (E2), the equation reduces to a linear equation:
n+m

n+m

j =1

j =1

lnαi + ∑ g ij ln X j = lnβi + ∑ hij ln X j

,

or
n

∑ ( g ij − hij ) ln X j = ln
j =1

Substituting ln X j = y j ,

ln

n+m
βi
− ∑ ( g ij − hij ) ln X jαij = n +1(E7).

βi
= bi . The equations (E7) is changed to the vector-matrix expression:
αi

r
r
r
(G D − H D )y D = b − (G I − H I )y I

(E8).

We substitute the matrixes as follows,

G D − HD = AD ,G I − HI = AI(E9)

to obtain the equation:rrr
AD y D = b − AI y I(E10),

where AD is the coefficient matrix of the kinetic order regarding dependent variables, AI is the coefficient matrix ofrr

the kinetic order regarding independent variables, y D is the solution vector of dependent variables, y I is thervector of independent variables, and b is the vector regarding the rate constants. Consequently, the vector of thedependent variables is symbolically solved as:rrr
y D = AD−1b − AD−1 AI y I(E11).

2.2 Sensitivity analysis
We define five types of gains and sensitivities. Here we explain the mathematical definition of them.
2.2.1. Logarithmic gain of a metabolite
In order to answer the question of how a relative change in an independent variable affect the steady-staterr

concentration of a metabolite, the logarithmic gain of a metabolite L( X D , X I ) is defined in terms of thecoefficient matrixes A −D1 and A I as:r r
L( X D , X I ) = − A −D1 A I(E12).
44
KIT, Kurata Laboratory2.2.2. Logarithmic gain of a flux
In order to answer the question of how a relative change in a flux affect the steady-state concentration of a

r

r

−1

metabolite, the logarithmic gain of a flux L(VD , X I ) is defined in terms of the coefficient matrices A D

and A I as:r rr r
L(V + , X I ) = G I + G D L( X D , X I )(E13).

2.2.3. Sensitivity with respect to a rate constantIn order to answer the question of how a relative change in a rate constant affect the steady state concentration of ametabolite, we solve the sensitivity with respect to rate constants, which is defined as:∂ ln X i∂ ln α j(E14),∂ ln X i∂ ln β j(E15).

S( X i , α j ) = −andS( X i , β j ) = −These sensitivities are expressed using vectors and matrixes as:r r
S ( X D , α) = − AD−1

(E16),r r
S ( X D ,β) = AD−1

(E17).

2.2.4. Sensitivity of dependent variables with respect to a kinetic order
In order to answer the question of how a relative change in a kinetic order affect the steady-state concentration of ametabolite, we define the sensitivity with respect to the kinetic order as follows,

S( X i , g ij ) = −∂ ln X i∂ ln g ij(E18),S( X i , h ij ) = −∂ ln X i∂ ln h ij(E19).

The partial derivatives of Eq. (E8) regarding gij are provided by:

or

r
∂AD r
∂y D
y D + AD
= 0,
∂g ij
∂g ij

(E20)

r
∂y D
∂A r
= − AD−1 D y D .
∂g ij
∂g ij

(E21).

Therefore,
45
KIT, Kurata Laboratory

r
r
∂A r
∂y D
= −g ij AD−1 D y D
S ( X D , g ij ) = g ij
∂g ij
∂g ij
r
r
∂A r
∂y
S ( X D , h ij ) = h ij D = −h ij AD−1 D y D
∂h ij
∂h ij

(E22)
(E23).

2.2.5. Sensitivity of a flux to a rate constantIn order to answer the question of how a relative change in a rate constant affect the steady-state concentration of aflux, we define the sensitivity with respect to a rate constant. We take logarithms in the first term of the right handside of S-system:n+mln Vl = lnαl + ∑ g lk ln X k

(E24)

k =1

Here, we differentiate to both side regarding lnαi to obtain:n

S (Vl ,αi ) =δil + ∑ g lk S ( X k ,αi )

(E25)

k =1

In analogy, partial derivatives of Eq.(E25) with respect to lnβi give:n

S (Vl ,βi ) = ∑ g lk S ( X k ,βi )(E26).

k =1

2.2.6. Sensitivity of a flux with respect to a kinetic order
In order to answer the question of how a relative change in a kinetic order affect the steady-state concentration of aflux, we define the sensitivity with respect to a kinetic order. The partial derivatives of Eq. (E24) with respect of ln
gij are provided as:n

S (Vl , g ij ) =δil g ij ln X j + ∑ g lk S ( X k , g ij )(E27).

k =1

In analogy, the partial derivatives of Eq. (E24) with respect of ln hij are provided as:n

S (Vl , h ij ) =δil hij ln X j + ∑ hlk S ( X k , hij )(E28).

k =1

2.3 Stability analysis
Defining the equation:
δXk = Xk － XkS (XkS steady state),
we apply the first order approximation to Eq. (E2) to obtain:n+m

g ikαi ∏ X j ijn

δX& i = ∑
k =1

j =1

Xk

n+ m

hikβi ∏ X j ijng

δX k − ∑X k = X kS

k =1

j =1

Xk

h

δX k (E29),

X k = X kS

orng ikViS+h V−
δX k − ∑ ik iSδX k
k =1 X kS
k =1 X kSnδX& i = ∑(E30),
46
KIT, Kurata Laboratory

As ViS+ = ViS− , we define:

V iS+
V iS−
=
X kS
X kS

f ik =(E31)

to deduce to:nδX& i = ∑ ( gik − hik ) f ikδX k(E32).

k =1

Stability at the steady state is investigated by solving the eigenvalues of the coefficient matrix consisting of:( gik − hik ) f ik

(E33).

3 Process flow
Figure D1 shows the whole processes of the simulator, where the process flow of S-system mainly consists of threeparts:
1) creation of the S-system coefficient file, which consists of the kinetic orders and rate constants that are solved insymbolic form at the steady state,2) creation of the S-system parameter file (the coefficients and steady-state concentrations) necessary forsensitivity/stability analysis, which is made using the S-system coefficient function and the steady state
concentrations.
3) simulation of S-system.
The checkdae module carries out the process of 1), following that the user function for ordinary differential
equations (TT, CMA, GMA, MM) are generated. However, the checkdae module does not convert the differential
and algebraic equations (DAEs), because they cannot be converted into S-system directly. The checkdae module
carries out the analysis of the sensitivities and stability of S-system by linking the S-system coefficient function to
the steady-state concentrations. The concentrations at the steady state must be solved before this analysis.

47
KIT, Kurata Laboratoryparsedae

DAEs(checkdae file)

Editionprocess

Edited DAEs
Filecheckdae

User Fucntion

Parameters

Sensitivity/stability analysis
S-system paremters

SolverDynaics

S-system coefficient
function

Steady state
anaylisSensitivity/stabil
ity results

S-systemparameters

S-systemuser function

solverS-systemdynamics/steady state
analysisFig. E1 Process flow for S-system analysis

48
KIT, Kurata Laboratory

4 Function
4.1 Creation of S-system parameter file
Following calculating the coefficients of S-system, the S-system parameter file is created. The S-system parametersprovide the values of the coefficients at the steady state for the S-system user function.

4.2 Sensitivity analysis
The coefficient matrix of S-system is created to solve the various kinds of the sensitivities.
4.3 Stability analysis
After the coefficient matrix of the system is created, the eigenvalues of the matrix are calculated by the Hessenberg
QR method.
5 Input/Output Specification

5.1 Input specification
The input files necessary for calculating the sensitivity and stability are the S-system user function that has beenobtained at the steady state, and the S-system parameter file. The S-system parameter file is as follows,##### Model ################
N_VAR;
3; # num of variables(all)
N_ALGEBR;
0; # num of variables(all)#Y_START ; index; initial_value; tag #comment
Y_START ;1; 6.7000e-02; X1#
Y_START ;2; 4.6501e-01; X2#
Y_START ;3; 1.5000e-01; X3##PARAM;name;index;val;val_start;num_survey;D/R/S;change-val;GA_min;GA_max;tag #comment
PARAM;constantPlayer;
1; 1.0000e+01;;0;D;;;; X4#PARAM;constantPlayer;
2; 5.0000e+00;;0;D;;;; X5#PARAM;constantPlayer;
3; 3.0000e+00;;0;D;;;; X6#PARAM;constantPlayer;
4; 4.0000e+01;;0;D;;;; X7#PARAM;constantPlayer;
5; 1.3600e+02;;0;D;;;; X8#PARAM;constantPlayer;
6; 2.8600e+00;;0;D;;;; X9#PARAM;constantPlayer;
7; 4.0000e+00;;0;D;;;; X10#
PARAM;kplus;1; 7.7884e-02;;0;D;;;; alpha_1#
PARAM;kplus;2; 5.8501e-01;;0;D;;;; alpha_2#
PARAM;kplus;3; 7.9346e-04;;0;D;;;; alpha_3#
PARAM;kminus;1; 1.0627e+00;;0;D;;;; beta_1#
PARAM;kminus;2; 7.9346e-04;;0;D;;;; beta_2#
PARAM;kminus;3; 1.0588e+00;;0;D;;;; beta_3
#
PARAM
;g
;
1; 6.6000e-01;;0;D;;;; gi_1_1
#
PARAM
;g
;
2; 1.0000e+00;;0;D;;;; gi_1_3
#
PARAM
;g
;
3; 9.5000e-01;;0;D;;;; gd_2_1#---omission--PARAM;h;5; -3.0600e+00;;0;D;;;; hd_2_3
#
PARAM
;h
;
6; 1.0000e+00;;0;D;;;; hi_2_5
#
PARAM
;h
;
7; 3.0000e-01;;0;D;;;; hd_3_3
#
PARAM
;h
;
8; 1.0000e+00;;0;D;;;; hi_3_6##T_EVENT ; name; index; time; value; #comment######## Don't edit the following lines. ######PARAM_NAME ;constantPlayer #PARAM_NAME ;kplus#
PARAM_NAME ;kminus#

49
KIT, Kurata Laboratory

PARAM_NAME ;g#PARAM_NAME ;h#5.2 Output specification
Following the analytical solution, the various sensitivities and stability are provided as follows.
##### Fri Mar 28 16:14:40 2003 ###### Logarithmic Gains and Sensitivities.
# The following data is the things when calculating STD.
##### Fri Mar 28 16:14:39 2003 #####
Solver No.
:
11
Param Survey:0
NR_TRIAL:20
NR_TOL_F: 1.000E-15
NR_TOL_X: 1.000E-15---omission--h[ 8]= 1.0000e+00#hi_3_6# Below this line is as a result of sensitivity analysis etc.
$ Fluxes $index,
y value ,1,
6.700e-02,2,4.650e-01,
3,1.500e-01,

Flux
1.068e+00, #X1
1.714e+00, #X2
1.714e+00, #X3$ Eigenvalues $No. ,Real ,
1, -2.508e+01,
2, -1.037e+00,
3, -5.281e+01,

Imaginary
0.000e+00
0.000e+00
0.000e+00$ Logarythmic Gains of Metabolites $param/var,y[ 1],
y[ 2],
y[ 3]constantPlayer[ 1], 8.283e-01, 1.029e+00, 1.216e+00
constantPlayer[ 2], 3.100e-01, 8.038e-01, 9.497e-01
constantPlayer[ 3], 1.255e+00, 1.560e+00, 1.843e+00
constantPlayer[ 4], -6.545e-01, -2.298e-03, -2.716e-03
constantPlayer[ 5], -8.648e-02, -2.243e-01, 3.264e-02
constantPlayer[ 6], -8.821e-01, -2.288e+00, -3.000e+00
constantPlayer[ 7], 3.681e-01, 9.545e-01, 1.128e+00$ Sensitivities of Metabolites with respect to rate constants $param/var,y[ 1],
y[ 2],
y[ 3]alpha[ 1], 1.255e+00, 1.560e+00, 1.843e+00
alpha[ 2], 9.686e-01, 2.512e+00, 2.968e+00
alpha[ 3], 8.821e-01, 2.288e+00, 3.000e+00$ Sensitivities of Metabolites with respect to kinetic orders $param/var,y[ 1],
y[ 2],y[ 3]g[ 1], 1.907e+00, 2.370e+00, 2.800e+00
g[ 2], 1.379e+00, 1.713e+00, 2.024e+00- - - - - omission - - - - h[ 7], 5.020e-01, 1.302e+00, 1.708e+00
h[ 8], -9.269e-01, -2.404e+00, -3.153e+00$ Logarythmic Gains of Fluxes $param/flux,
v[ 1],

v[

2],

v[

3]

50
KIT, Kurata LaboratoryconstantPlayer[
constantPlayer[
constantPlayer[
constantPlayer[
constantPlayer[
constantPlayer[
constantPlayer[

1],
2],
3],
4],
5],
6],
7],

6.600e-01, 3.649e-01, 3.649e-01
0.000e+00, 2.849e-01, 2.849e-01
1.000e+00, 5.528e-01, 5.528e-01
0.000e+00, -8.147e-04, -8.147e-04
0.000e+00, 9.791e-03, 9.791e-03
0.000e+00, 9.987e-02, 9.987e-02
0.000e+00, 3.383e-01, 3.383e-01$ Sensitivities of Fluxes with respect to rate constants $param/var,
v[ 1],
v[ 2],v[ 3]alpha[ 1], 1.000e+00, 5.528e-01, 5.528e-01
alpha[ 2], 0.000e+00, 8.903e-01, 8.903e-01
alpha[ 3], 0.000e+00, -9.987e-02, 9.001e-01$ Sensitivities of Fluxes with respect to kinetic orders $param/var,
v[ 1],
v[ 2],
v[ 3]g[ 1], 1.520e+00, 8.401e-01, 8.401e-01
g[ 2], 1.099e+00, 6.073e-01, 6.073e-01
g[ 3], 0.000e+00, -2.286e+00, -2.286e+00
g[ 4], 0.000e+00, 2.795e-01, 2.795e-01- - - - - omission - - - - h[ 6], -3.170e-16, 5.387e-01, -4.374e+00
h[ 7], -1.545e-17, -5.684e-02, -5.684e-02
h[ 8], 2.851e-17, 1.049e-01, 1.049e-01

51
KIT, Kurata Laboratory

F. OPTIMIZER (GA)
1. Introduction
The Genetic Algorithms (GAs) are known as one of the algorithms that can seek out the global minimum. GA isbased on the heuristic assumptions that the best solution will be found in the regions of the parameter spacecontaining a relatively high proportion of good solutions and that these regions can be explored by the genetic
operators of selection, crossover, and mutation. GA offers a number of advantages, e.g., GA searches form a set ofdesigns, and can explore the parameter space without trapping in local optima. However, this method has the
disadvantage that the computational cost of many runs for the design code is considerably large, and there havebeen many parameters regarding the operations to determine with experiences and intuitions. In order to execute
GAs efficiently and systematically, we have developed the search modules The algorithms suitable for optimizing
biochemical networks are presented so that users can make a efficient strategy for the parameter search. This search
can be readily carried out on parallel computers.

2. Function specification and application problems
2.1 Function Specification
Algorithm- Encode method-Real-coded GA (RGA)-Bit-string GA (BGA)-Binary coding-Gray coding- GA type[RGA, BGA]-Distributed GA-Distributed and Integrated GA-Generation alternation[RGA, BGA]

-Elite conservation strategy (Normal)-Minimal Generation Gap (MGG)-Selection method:

Roulette, Tournament, Random.

- Crossover-BLX (Blend Crossover)[RGA]-UNDX (Unimodal Normal Distribution crossover)[RGA]-UNDXm (Multi-parental Unimodal Normal Distribution crossover) [RGA]-SPX (Symplex crossover)[RGA]-n-point crossover [BGA]- Mutation-Uniform mutation within region [RGA]
52
KIT, Kurata Laboratory

-Uniform mutation with fixed width [RGA]-Normal mutation with fixed width [RGA]-Uniform mutation with variable width [RGA]-Normal mutation with variable width [RGA]-Bit reverse mutation [BGA]
Input files and data for GAs- GA-parameter setting file
- GA-population setting file- parameter file (created by the checkdae module)- parameter-range setting file- Fitness function
Output file and data
Parameters for executing GAs (standard output, output file)
Information regarding the best individual of each generation(standard output)

Information regarding the whole individuals every generation (output file)

2.2 Application problem
This search engine module can be applied to a maximization problem that defines a fitness function by using
continuous and independent variables with the limited ranges.
Exceptionally, this search engine can be applied to the optimization of discrete problems such as a knapsackproblem. In the bit string GA for which the range of variables is set [0, 1] and the quantization number is 1, thevariables are determined either 0 or 1.
Note that the range is closed.

3. Execution of programs
Execution on a command line
% load module<arguments>

Arguments:-r[parameter-range setting file]or –I [parameter file] (required)-p[GA-parameter setting file](required)-s[GA-population setting file](option)-f[Output file]-hShow help.

The parameter file is exactly the same as the parameter file that the CADLIVE Simulator creates. If users employ
the parameter file to determine the regions of the parameters, the parameter-range setting file is not required.
53
KIT, Kurata Laboratory

Examples:(1)%load_module_name–r domain.dat -p param.dat -f output.dat -s init.dat(2)%load_module_name-h

For MPI,(3)%mpirun –np <number of process> <load module name><arguments>%mpirun –np 3 <load module name> –r domain.dat -p param.dat-f output.dat

For MPICH-G,(4)%mpirun –globusrsl <rsl file>The "rsl" file is described in the manual of MPICH-G.
One process for one CPU is recommended.
How to connect the GA module to the CADLIVE Simulator
Users are able to connect the GA module to the CADLIVE Simulator, where the ranges of the parameters aredefined in the parameter file (parameter setting file). In order to use the parameter file as the input file for GA
modules, users have to edit the parameter file whose format is suitable for the GA module. The program is executed
as follows:%load module name- I[parameter file]-p[GA parameter setting file]- f[GA result file].

The detailed instruction is described elsewhere.

4. Input/output specification
4.1 Parameter-range setting file
The closed ranges of the parameters are defined in the parameter-range setting file, and the number of the records isthe number of the parameters to search, as shown in the following box.
# Format for a parameter range-setting filelower value of range, upper value of range

To maximum of 20 records.
# comments

Format: Fields are sectioned by using comma (,). The preceded parts by # are comments.
54
KIT, Kurata Laboratory

For example, the ranges of three parameters are [-5.12,5.12], [-5.12, 0.0], [0.0, 3.14] as follows.
# Example: parameter range-setting file
-5.12, 5.12

# x1

-5.12, 0.0# x20.0, 3.14# x3

4.2 GA-parameter setting file
4.2.1 GA-parameter setting file
The control parameters for GAs are classified by five key words: encode method, GA type, digenesis (generation
alternation), crossover, and mutation. The GA-parameter setting file is provided as follows.
# Format of GA-parameter setting file
Key word; Alternatives; parameter_1; ****; parameter_n;# comments

The preceded parts by # are comments. The key words and alternatives are selected from the strings shown in Table
F1.
Table F1 Control parameters for setting GAs. Key words and their alternatives are selected. Upper and lower cases
are distinguished for defining the parameters. The program runs until the fitness becomes more than the value forterminating a search. Both RGA and BGA can select the transparent alternatives. The thin gray alternatives can beselected only by RGA, and the dark gray alternatives can be selected only by BGA.
Key word
ENCODE(Encode method)

Alternatives
REAL
BINARYGRAYGATYPEDGA(Island model)
DIGA
DIGENESISNORMAL(Generation alternation) MGG
CROSSOVER
BLX(Crossover method)
UNDX
UNDXm
SPX
NPOINTSNONE
MUTATION
RegionUni(Mutation method)
FixedUni
FixedNormal
VariableUni

Meaning
Real GA
Binary coding bit string type GA
Gray coding bit string type GA
Distributed GA (island model)
Distributed and Integrated GA
Normal generation
MGG
Blend crossover
Unimodal Normal Distribution crossover
Multi-parental Unimodal Normal Distribution crossover
Symplex crossovern-point crossover (BGA)
No crossover
Uniform mutation within region
Uniform mutation with fixed width
Normal mutation with fixed width
Uniform mutation with variable width
55
KIT, Kurata Laboratory

VariableNormal
BitReverse
NONE

Normal mutation with variable width
Bit reverse mutation
None

Users can set the parameters with respect to each alternative in the following box.
# Parameter description with respect to key words and alternatives.
ENCODE; REAL;ENCODE; BINARY; quantization number;
ENCODE; GRAY; quantization number;
GATYPE; DGA; maximum generation number; value for terminating a search; number of islands; population
number within islands; immigration method; immigration interval; immigration rate;(The final two parameters are required for immigration operation. Actually the above parameters must be on one
line.)GATYPE; DIGA; maximum generation number; value for terminating a search; number of islands; population
number within islands; generation number for integration;
DIGENESIS; NORMAL; number of elites; rule for selection; tournament size (if tournament is employed);
DIGENESIS; MGG; number of children generated;
CROSSOVER; BLX; alpha;
CROSSOVER; UNDX; alpha; beta;
CROSSOVER; UNDXm; alpha; beta; m;
CROSSOVER; SPX; epsilon;
CROSSOCER; NPOINTS; n;
CROSSOVER; NONE;
MUTATION; FixedUni; mutation rate; width of parameter range;
MUTATION; FixedNormal; mutation rate; standard deviation;
MUTATION; VariableUni; mutation rate;
MUTATION; VariableNormal; mutation rate;
MUTATION; RegionUni; mutation rate;
MUTATION; BitReverse; mutation rate;
MUTATION; NONE;
In the bit-string GA, the real values are converted into the bit sequences by dividing the range of the variables with

2 n . The number of n that determines the resolution is named as the quantization number. A part of the parameters
for GA type, generation alternation, crossover, and mutation depends on the parameters of the encode methods. The
parameters for GA type, generation alternation, crossover, and mutation are independent mutually.

56
KIT, Kurata Laboratory

Table F2 Values that users are allowed to set with respect to each parameter. Note that the GA module candistinguish the upper or lower cases of strings.
Parameter
quantization number
maximum generation number
value for terminating a searchnumber of islands
population number within
islands
immigration switch
immigration interval
immigration rate
generation for integration
number of elites
selection rule

ValuesDefault value
1
Integer ≥ 1
100
Integer ≥ 1
real value (double type)-1.00E-10
3Integer ≥ 1variable number +2 ≤ Integer
30≤ Maximum
ON, OFFON
3
Integer ≥ 1
Real value [0,1]15Integer ≥ 1
Integer ≥ 0
population number
Roulette, Tournament, Random
Tournament
1 ≤ Integer ≤ (population number
2
size of tournamentof islands – number of parents +1)
number of children generated
population
1 ≤ Integer ≤ Maximum
by MGG
number of islands
alpha (BLX)
Real value > 00.5
alpha (UNDX)Real value > 0
1alpha (UNDXm)
Real value > 0
1beta (UNDX)Real value > 0
0.35
beta (UNDXm)
Real value > 00.35
max(1, the number of variables
1 ≤ Integer ≤ number ofM-2)variables
1/sqrt (number of variables
epsilon
Real value > 0+2)N1
Integer ≥ 1
mutation rate
Real value [0,1]
1/ number of variables
widthReal value > 0
20
standard deviation
Real value > 02

4.2.2 Examples for GA parameter setting file
Example 1:
Real GA
Distributed GA, To the maximum of 100, Value for terminating a search = 0, number of islands = 5, population
number within islands = 100, Immigration is ON, Immigration internal = 10, Immigration rate = 0.25, MGG,
children = 10, UNDX, parameter defaults, No mutation is applied.
# example 1
ENCODE; REAL;GATYPE; DGA; 100; 0; 5; 100; ON; 10; 0.25;
DIGENESIS; MGG; 10;
57
KIT, Kurata Laboratory

CROSSOVER; UNDX; ; ;MUTATION; NONE;
Example 2
Binary coding bit-string GA, Quantization number = 10, Distributed and integrated GA, maximum generation 200,
value for completing execution = 0, number of island = 3, population number within island = 30, generation number
for integration = 20, Normal generation alternation, number of elites = 5, Roulette selection, Two-point crossover,
Bit-reverse mutation, mutation rate 0.3# example 2
ENCODE; BINARY; 10;
GATYPE; DIGA; 200; 0; 3; 30; 20;
DIGENESIS; Normal; 5; Roulette;
CROSSOVER; NPOINTS; 2;
MUTATION; BitReverse; 0.3;
4.2.3 GA-population setting file
The GA module is able to read the result file from the random search. In the n-dimensional parameter space to be
searched, the GA-population setting file is provided by:# Format for GA-population setting file
SearchID, fitness, parameter1, ....., , parameter_n....................................
**
SearchID, fitness, parameter1, ....., , parameter_n
The symbol of "#" indicates that the line is comments. The section between islands is expressed by "**". Users
neither have to put "**" on the initial island nor under the final island. If such descriptions are employed, an empty
island is regarded to exist.
4.2.4 Examples for GA-population setting fileThree islands have five individuals, respectively# island No1
No0-67,-0.304018,0.915873,3.28519
No0-66,-9.88925,2.14165,1.47417
No0-43,-21.9843,2.8856,1.10932
No0-51,-25.6788,2.53911,1.00846
58
KIT, Kurata Laboratory

No0-5,-28.4887,0.792621,2.64166**#island No2

No0-53,-3.04234,0.987835,2.93188
No0-41,-17.3451,1.64468,1.63063
No0-71,-23.0538,1.26285,1.99734
No0-40,-25.7664,1.81623,1.37769
No0-42,-29.2198,2.79655,1.15365**#island No3

No0-10,-4.99208,0.972839,3.17752
No0-88,-17.9748,1.64134,1.6242
No0-44,-24.7664,1.06062,2.22619
No0-16,-27.1525,1.76528,1.33739
No0-6,-30.4514,1.33549,1.56515

4.3 Standard output
Before execution, the values of the parameters are displayed on the screen as follows.
===GA PARAMETERS======================
Encode Type: Real-coded GA

GA type: DGA

Max Generation: 1000Fin value: -1.000000e-10

Total population: 1000

Number of islands : 10
Island population: 100

Digenesis type
Number of elites: NORMAL:1

Selection type: Tournament

Tournament Size

:2

Crossover type

: SPX

SPX(epsilon): 2.645751

Mutation type

: RegionUni

Mutation rate

: 0.200000
59
KIT, Kurata Laboratory

Immigration type

: Random ring

Immigration interval:5

Immigration rate

: 0.500000===Searching space information==========
Number of Variable:5

Domain variable[0]: [-5.120000, 5.120000]

Domain variable[1]: [-5.120000, 5.120000]

Domain variable[2]: [-5.120000, 5.120000]

Domain variable[3]: [-5.120000, 5.120000]

Domain variable[4]: [-5.120000, 5.120000]========================================A summary for the GA-population setting file is displayed, which includes the file name, the island numbers, thepopulation number within each island, the number of short islands, and the number of redundant islands. In the caseof the short islands, the program generates a new island with a random population automatically. When the number
of islands is larger than the appointed one, the number of the redundant islands is displayed.
===Initial Generation File==============
File name:rs_out_list
Island, Population, # of read, Random, Redundant
Island[0], 15, 20, 0, 5Island[1], 15, 19, 0, 4
Island[2], 15, 5, 10, 0
Island[3], 15, 31, 0, 16
Island[4], 15, 11, 4, 0
Number of redundant islands: 2#

displayed when redundant islands exist.

========================================The information regarding the best individual is displayed every generation on the screen. After completing thesearch, the result is displayed. The bit-string GA indicates the real value sequences, not the bit sequences.
Indication of results every generation:
generation[0], Obtain: 0, Fitness: -26.652798, Gene: (-1.982981, -4.013446, 0.087129, -1.089640, 0.110301)
generation[1], Obtain: 1, Fitness: -22.533149, Gene: (-0.890659, -1.025080, -2.211625, 0.958039, -1.936774)
generation[2], Obtain: 2, Fitness: -14.719590, Gene: (0.074704, -0.078463, -0.083288, 1.033160, -2.918250)(omission)(when the target fitness or the maximum generation number is achieved.)******Search Result******
60
KIT, Kurata Laboratory

Best geneObtained Generation: 65

Fitness: -0.000000

Gene: (-0.000000, -0.000000, 0.000000, 0.000000, 0.000000)*************************

4.4 Output file
The output file is appointed as an argument on a command line. The results regarding the entire individuals every
generation are written in the form of "csv".
generation,0Best gene:Obtained generation, fitness, variable[0], variable[1], variable[2], variable[3], variable[4]0, -46.007195, -3.995016, -0.149753, 3.052604, 3.074872, -2.078418
ID, fitness, variable[0], variable[1], variable[2], variable[3], variable[4][0][0], -46.007195, -3.995016, -0.149753, 3.052604, 3.074872, -2.078418[0][1], -65.409291, -2.949895, 4.444091, 1.857122, 1.830248, -1.036918[0][2], -83.227264, 1.431177, 2.992203, 3.875130, 4.705202, 0.037506(omission)[1][7], -100.625396, -4.595544, 3.022847, -4.746842, -3.206425, -0.934593[1][8], -110.044951, 2.700800, 4.610274, -0.096368, -4.767272, 1.675321[1][9], -135.758812, 3.813671, -4.647205, -3.593073, -4.125371, 5.059261
The brackets [][] indicate [island number][individual number within islands]. The individuals within each island aresorted according to the values of the fitness. A smaller number indicates a higher fitness.

5. Detailed specification of functions
5.1 Encode method
This optimization program presents the real-coded GA (RGA) and the bit-string GA (BGA) as the method that
encodes design variables. For the bit-string GAs, the binary coding or gray coding is employed as the method thatexchanges between integers and bit-string (0,1).
5.1.1 Bit-string GA
The bit-string GA encodes the design variables as the bit string that consists of 0 and 1. The bit-string GA divides
the range of variables into 2

nregions, and assigns the integer of 0, K ,2 − 1 for each region in a numericalnorder, thereby converting the integer to the bit-string. Two encoding methods as follows:
Binary coding: converting an integer by binary arithmetic.
Gray coding: converting an integer with a binary code under the condition that the Haming distance between
two codes to continue is always kept 1. The Haming distance is defined as the number of the figures with
61
KIT, Kurata Laboratorydifferent values.
Table F3

Integer
0
1
2
3
4
5
6
7

Binary and gray codes

Binary code
000
001
010
011
100
101
110
111

Gray code
000
001
010
011
110
111
101
100

5.1.2 Real GA
The real GAs employ the sequence of real values as the chromosome.

5.2 GA typeThis module presents two types of GAs, the distributed GA (DGA, island model) and the distributed and integrated
GA (DIGA).
5.2.1 Distributed GA
The distributed GA divides the population into multiple islands. The individuals of each island evolve within theisland, and some individuals immigrate mutually among the islands. The immigration method can be selected fromrandom swap, random ring, and no immigration.
Random swap: selecting immigrants randomly out of randomly selected islands under the given immigration
rate in order to exchange two individuals between the islands. The elites of each island are not allowed toimmigrate.
Random ring: determining the destination of immigrants so that the way between the destination and origin of
the immigrants forms a ring. The elites of each island are not allowed to immigrate.
5.2.2 Distributed and integrated GA (DIGA)The DIGA divides the population into multiple islands as well as DGA. Each island evolves independently until all
the islands are integrated to a new island at the generation number for integration. At their integration, the elites ofeach island are gathered to the new one, and the remaining individuals are selected randomly. The DIGA is notallowed to operate the immigration.

5.3 Generation alternation
The GA module presents two methods for generation alternation, the ordinary generation alternation method(employing elite conservation method), and Minimal Generation Gap (MGG). The ordinary method is sometimes
called the simple GA, but such a name may cause users to misunderstand the methods. Here, we name the ordinary
62
KIT, Kurata Laboratorygeneration alternation method.
5.3.1 Ordinary generation alternation method
The ordinary method replaces all the parents by children. Actually, the parents except the elites are replaced,because the elite conservation strategy is employed. If users want to replace all the parents including the elites, set“the number of elites” as zero.
The roulette, tournament, and random selections are presented as the methods for selecting the parents to crossover.
The roulette is made based on the difference:(difference)=(fitness of an individual) – (the lowest fitness of the island).

5.3.2 Minimal Generation Gap (MGG)
The MGG is the most desirable model that can avoid early convergence and suppress evolutionary stagnation. The
procedure of MGG can be summarized in the following box.
1. Create an initial population randomly
2. Sample two individuals randomly without replacement from the population
3. Generate children from the selected parents and characterize them
4. Replace two individuals (parents) by the best individual and the roulette-selected individual, resulting in the next
generation.
The above method is the original MGG. Some crossover methods require more than two parents. Thus, thecrossover is defined by the rules:
BLX-α, n -point mutation: The parents to replace are the same ones to crossover.
UNDX, UNDX-m, SPX: Two parents to replace are selected randomly out of all the parents that have been
employed to generate children.

5.4 Crossover method
The crossover method to apply depends on the encoding methods, the bit-string GA (BGA) and the real GA (RGA).
In this section, [RGA] indicates that the crossover is available for RGA, [BGA] for BGA.
5.4.1 BLX-α [RGA]The BLX is a simple algorithm among the real-coded GAs, but its performance is not guaranteed in thevariable-dependent functions.

ci = u (min( pi1 , pi2 ) − αI i , max( pi1 , pi2 ) + αI )I i =| pi1 − pi2 |,where p = ( p1 , K , p n ), p = ( p1 , K , p n ) are the parents, c = (c1 , K , c n ) is the child, and
1

1

1

2

2

2

u (⋅)

represents the function that generates a uniform random number within the closed region of the parentheses.
63
KIT, Kurata Laboratory5.4.2 UNDX / UNDX-m [RGA]The UNDX-m has been proposed to improve a search performance on an optimization problem with poorly scale
coordinate systems. The UNDX-m generates offspring vector values by sampling values from the m-dimensional
space that the m+1 parents span around their middle point. The UNDX is the same as the UNDX-m at m = 1 except
some respects. The prototype algorithm of the UNDX-m in the n-dimensional parameter space is provided by:
1(1) Select m +1 parents x , K , x

m +1

randomly from the population.

(2) Let the center of mass of these parents

1 m +1 ip=∑x ,m + 1 i =1
and let the difference vector between(3) Select another parent x

m+ 2

xi and p be d i = x i − p .

from the population randomly= x m + 2 − p orthogonal to d 1 ,K , d m .
1
n− m
1
m
be orthogonal bases of the subspace orthogonal to the vectors d , K , d .
(5) Let e , K , e(4) Let D be the length of the component of d

m+ 2

c

(6) Generate the child x by:mn−m

x = p + ∑ wi d + ∑ vi De i ,
c

i

i =1

i =1

where wi , vi are the random numbers that follow normal distributions N (0, σ ξ ), N (0, σ η ) , respectively. The
2

2

parameters, σξ and ση, are provided by:σ ξ = α / m,

ση =

1

m +1

3

n−m

m+2

2

β

where α= 1.0 and β= 0.35 are recommended.
The UNDX is carried out as follows:12(1) Select three parents x , x , x

3

randomly.

(2) Let the center of mass of these parents:x p = ( x1 + x 2 ) / 2 ,and let the difference vector d = x − x .
123(3) Let D to be the distance between the third parent x and the line connecting x1 to x2.
(4) Generate a child xc by the following equation:n −1

x = x + ξd + D ∑η i e i ,
c

p

i =1

where

ξ and η i are the random numbers that follow normal distributions N (0, σ ξ2 ) and N (0, σ η2 ) ,respectively. The parameters, σξ and ση, are given by:σξ = α,ση = β / n64KIT, Kurata Laboratorywhereα = 0.5 and β= 0.35 are recommended.
Difference between UNDX and UNDX-m
The recommended value of α is a half of that of the UNDX-m, because the UNDX employs the difference vector
between the parents instead of the difference vector based on the center of mass. Actually, both are the same. On theother hand, for the sub-search component, the UNDX-m provides:σ η = β / n − 1 at m = 1,
where the difference is indicated clearly between them.
Table F4 Difference between UNDX and UNDX-m
Revision on sub-search componentσξ = αUNDXUNDX-mση = β / nσξ = α / m

ση =

1

m +1

3

n−m

m+2

2

β (m = 1, β / n − 1)Note: Rank cancellation
In the UNDX-m, the space that the difference vectors span is the main search component, and the orthogonal
subspaces are the minor search components. If the difference vector sets that the main search component belongs tohave any dependent element, the dimension of the main search space is less than the number of the vectors thatspan it, canceling the rank.
When the rank is cancelled, the optimization module does not change the main search component, and lets theminor search component to be the subspace orthogonal to the space (P) that independent vectors of the mains
search components span. The parameter D is set as the component of the vector dm+ 2

orthogonal to the space (P).

5.4.3 SPX [RGA]
SPX generates offspring vector values by uniformly sampling values from the simplex formed by multiple parent
vectors.

rr(1) Seletct m+1 parents P0 , K , Pn randomly from the population.

r(2) Let the center of mass of the parents G ben rr
G = ∑ Pi .
i =0

r

r

(3) The vector x k , C k are determined as follows:

65
KIT, Kurata Laboratory

r
r
r
r
x k = G + ε ( Pk − G ) (k = 0, K, n)rr,
0(k = 0)r
Ck = 
rr
rk −1 ( x k −1 − x k + C k −1 ) (k = 1, K n)
where ε is the extension coefficient (ε > 0), rk is provided by random uniform function:

rk = (u (0,1))1 /( k +1)(k = 0, K, n − 1) .

The coefficient ε is provided by:

ε = n+2(4) The child is provided byr rr
C = xn + Cn

5.4.4 N point crossover [BGA]
N points to crossover are selected at a chromosome randomly, and the parameters between the crossing points are
replaced alternately by the corresponding parameters of a parental chromosome.

5.5 Mutation method
The mutation methods to apply depend on BGA or RGA. In this section, the mutation method for RGA is explained.
The individuals that have been selected out of an island at a given mutation rate are mutated according to thevarious methods, where ( x1 , K , x n ) is the coordinate of an individual, and ( z1 , K , z n ) is the coordinate thatshows the best individual in all of the individuals within the island.
5.5.1 Uniform mutation within the region [RGA]
The vector of the selected individual is changed under the uniform random distribution within the given region.
5.5.2 Uniform mutation with fixed width [RGA]
A mutated vector is given by the function:u ( xi − w, xi + w) ,which generates a uniform random number between xi − w and xi + w , where w is the width that users givearbitrarily.
5.5.3 Normal mutation with fixed width [RGA]
A mutated vector is provided by the function:u ( xi − w, xi + w) ,which generates a normal random number between xi − w and xi + w , where w is the width that users givearbitrarily. The coordinate of the mutated vector follows a normal distribution:
66
KIT, Kurata LaboratoryN ( xi , s 2 ) ,where s is the standard deviation that users can give arbitrarily.
5.5.4 Uniform mutation with variable width [RGA]
A mutated vector is provided by the function:u ( xi − wv , xi + wv ) ,which generates a uniform random number between xi − wv and xi + wv , where wv is given by:

wv =| z i − xi | .
5.5.5 Normal mutation with variable width [RGA]
A mutated vector is provided by the function:u ( xi − w, xi + w) ,which generates a normal random number between xi − w and xi + w , where w is the width that users givearbitrarily. The coordinate of the generated vector follows a normal distribution:N ( xi , s v2 ) ,where the standard deviation is s v =| z i − xi | .
5.5.6 Bit reverse mutation [BGA]
A vector is mutated by reversing the bit according to the given mutation rate.

6. Flow chart
6.1 Single without MPI+-main+-initialize parameters+-display the employed parameters+-generate the initial population and characterize the fitness+-if(ORDINARY )+-NormalExecute()+-for (final generation)+-if (DGA && immigration generation) immigrate+-if (DIGA && integration generation) integrate+-for (all the islands)+-for (non-elite individuals)+-select parents based on the crossover method+-generate children by crossover+-replace parents by children (except elites)

67
KIT, Kurata Laboratory

+-for (non-elite individuals) mutate+-characterize the fitness of non-elites+-sort the individuals according to their fitness+-conserve the best individual of all the islands+-if (termination condition) break;+-display the results+-else if (MGG)+-MGGExecute()+-for(final generation)+-if (DGA && immigration generation) immigrate+-if (DIGA && integration generation) integrate+-for (all the islands)+-sampling parents without replacement+-generate a specific number of children+-characterize the fitness of children+-select the best individual and another by roulette selection+-replace parents by children+-conserve the best individual of all the islands+-if (termination condition) break;+-display the results

6.2 MPI
6.2.1 Master-slave modelThis optimization module presents a master-slave model for parallel computing, where the master manages thewhole processes and send jobs to the slaves. Concretely speaking, the master manages the whole population, anddisplays the results. On the other hand, the slaves execute the evolution within the islands.
Master:
Generating the whole population
Sending jobs to slaves.
The queue manages jobs and asks the slave that has terminated the task to do another job
Managing immigration or integration
Displaying the results
Slaves:Independently processing a job
Crossovering and mutating chromosomes
Characterizing the fitness and sorting the individuals.
68
KIT, Kurata Laboratory

The speed of parallel computation mainly depends on the grain of calculation provided for each CPU. The small
grain increases the overhead of communication, resulting in decreasing the advantage of parallel computation.
Delivering a job to a slave by an island unit increases the grain of calculation. The calculation rate of a simple
parallel system is limited by the slowest CPU. By contrast, the use of the master-slave model loads many jobs to the
CPU with high-performance, thus the calculation rate is enhanced more than the simple parallel system.
The master-salve model is shown in the following diagram.

Masterdelivering a job (Island) to a slave
operating immigration
output the results
Island1IslandIsland

3

1

2

Slave
crossovermuatatonfitness calculation
sort

IslandIsland

1. ask a job to slave

Slave

2. return to a signal

crossovermuatatonfitness calculationsort

3. ask a new job tothe salve that returned
a signal

6.2.2 A flow char in MPI+-main+-initialize MPI+-if (MASTER)+-input the GA parameter setting file and variable range setting file+-send the data on parameter and range to all the slaves+-display the employed parameters+-if (SLAVES)+-receive the data from the master+-register the data+-if (MASTER)+-generate the initial population of all the islands+-for (SLAVES)+-send the data of the initial population+-if (All the islands are sent) break;+-for (infinite loop)

69
KIT, Kurata Laboratory

+-receive the response from a slave+-if (there is an island to send data) send the island to the slave+-if (All the islands are received) break;+-for (All the SLAVES) send the command to exit the loop+-if (SLAVE)+-for (infinite loop)+-receive data from the master+-if (Exiting the loop is commanded) break;+-characterize the fitness of all the individuals+-send the data to the master+-if (ORDINARY)+-normalExecute()+-if (MASTER)+-for (final generation)+-if (DGA && immigration generation) immigrate+-if (DIGA && integration generation) integrate+-for (All the SLAVES)+-send the data of all the islands+-if (All the islands are sent) break;+-for (infinite loop)+-receive the response from a slave+-if (there is a island to send)
send the island to the slave+-if (all the islands are received) break;+-conserve the best of all the islands+-if (termination condition)+-send the command for exiting the loop, break;+-display the results+-if (SLAVE)+-for (infinite loop)+-receive the data from the master+-if (Exiting the loop is commanded) break;+-for (non-elite individuals)+-select parents based on the crossover method+-generate a child by crossover+-replace parents by children (except elites)+-for (non-elite individuals) mutate+-characterize the fitness of non-elites

70
KIT, Kurata Laboratory

+-sort the individuals according to the fitness+-send data to the master+-else if (MGG)/* MGG */+-MGGExecute()+-if (MASTER)+-for (final generation)+-if (DGA && immigration generation) immigrate+-if (DIGA && integration generation) integrate+-for (All the SLAVES)+-send the population data of all the islands+-if (All the islands are sent) break;+-for (infinite loop)+-receive the response from a slave+-if (there is an island to send) send the island to the slave+-if (All the islands are received) break;+-conserve the best individual of all the islands+-if (termination condition) break;+-send the command for termination to the slaves,
break;+-display the results+-if (SLAVE)+-for (infinite loop)+-receive the data from the master+-if (Termination is commanded) break;+-select parents randomly+-generate a specific number of children from parents+-characterize the fitness of children+-select the best individual and another by roulette selection+-replace parents by children+-if (Mutation)+-mutate+-characterize the fitness and sort+-send data to the master.

71
KIT, Kurata Laboratory

G. CLIENT-SERVER MODEL
The CADLIVE Simulator is a client-server model, supporting Internet Explorer. Users access a server through php,
and operate on Web pages. The system functions on LINUX (Red Hat Version 7.1), and employs Xerces2 Java

Parser for parsing the regulator-reaction equations in XML, postgreSQL for managing database, the C
programming language for calculating differential and algebraic equations and matrixes, and Message-Passing
Interface (MPI) for parallel computing.

72
KIT, Kurata Laboratory