An improved cellular automata model of enzyme kinetics based on genetic algorithm

Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎ Contents lists available at ScienceDirect Chemical Engineering Science journal homepage: www.elsevier...

Download PDF

1MB Sizes 5 Downloads 32 Views

Report

PDF Reader
Full Text

Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

An improved cellular automata model of enzyme kinetics based on genetic algorithm Saurajyoti Kar a, Kaustuv Nag b,1, Abhishek Dutta c, Denis Constales d, Tandra Pal b,n a

Department of Biotechnology, National Institute of Technology, Durgapur 713209, West Bengal, India Department of Computer Science & Engineering, National Institute of Technology, Durgapur 713209, West Bengal, India c Department of Chemical Engineering, Groep T – Leuven Engineering College, (Associatie KU Leuven), A. Vesaliusstraat 13, B-3000 Leuven, Belgium d Department of Mathematical Analysis, Ghent University, Galglaan 2, B-9000 Ghent, Belgium b

H I G H L I G H T S

Cellular automata are used to model enzymatic reaction kinetics on a 2D torus. Fine-tuning the probabilities of reactions and movement is challenging. A multi-objective genetic algorithm obtains the probabilities in feasible time. A validation study indicates the applicability of the cellular automata model.

art ic l e i nf o

a b s t r a c t

Article history: Received 28 April 2013 Received in revised form 30 June 2013 Accepted 8 August 2013

Based on a generalized Cellular Automata (CA) model which has the ﬂexibility of reaction interactions between all possible two-agent combinations on both reactant and product side, in previous studies by the authors, probabilistic reaction events were combined in a multiplicative order and an exploratory study of the variability analysis was performed. However it was found that the probabilistic parameter values are extremely sensitive to the ﬁnal model output. As such, NSGA-II based on a Multi-objective Genetic Algorithm (MOGA) is used in this study to obtain optimum sets of these probability parameter values within a feasible computational time. For each generation of the MOGA, the probability rules parameter values are updated (improved) based on the objective functions. Finally, a validation study is performed to indicate the applicability of the CA model to represent speciﬁc enzymatic reactions when coupled with the multi-objective optimization algorithm. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Cellular automata Multi-objective genetic algorithm NSGA-II Enzyme kinetics

1. Introduction Natural processes from population dynamics (Soarse-Filho et al., 2002; Lauf et al., 2012) to ﬂuid dynamics (Somers, 1993; Chopard and Masselot, 1999), from bioﬁlm growth (Pizarro et al., 2005; Fagerlind et al., 2012) to tumor growth (Kansal et al., 2002; Alemani et al., 2012) all occur in discrete space and time so that a probabilistic nature of representation is needed to understand precisely the dynamics of such processes. Stochasticity is an intrinsic property of any natural system, more speciﬁcally at a microscopic scale (Deutsch and Dormann, 2005). Several studies have observed occurrence of stochastic ﬂuctuations in living systems (Abkowitz et al., 1996; Meng et al., 2004; Wu and Higgs 2012). Most biological events are triggered as network of chemical reactions which pose larger challenges due to computational n

Corresponding author. Tel.: þ 91 9434537021; fax: þ 91 343 2547375. E-mail address: [email protected] (T. Pal). 1 Present address: Instrumentation and Electronics Engineering, Jadavpur University, Kolkata 700 032, West Bengal, India.

representation and requirements (Kier et al., 2005). The noise of one event can alter the outcomes of the succeeding event and simultaneously affect the ﬁnal outcome (Blake et al., 2003). When concentrations of species in an enzymatic reaction are high, the reaction can be mathematically modeled using deterministic approaches, like Michaelis–Menten kinetics. But, at lower concentrations, it is advisable to use modeling tools which have some stochastic property within the rule structure. Such models shall evolve to show certain global properties and are also unpredictable until the event actually occurs to completion. Thus an accurate mathematical model can help clarify the roles of individual components within a process and generate speciﬁc, testable hypotheses and predictions (Apte et al., 2008). Cellular Automata (CA) is one such stochastic approach in which a natural event can be successfully simulated. Using probabilistic rule application, CA can capture the properties that evolve from the stochastic nature present in microscopic scale in natural systems assuming both time and space to be discrete. In our earlier studies (Kar et al., 2010; Dutta et al., 2011), a generalized two-dimensional (2D) interaction model was created to represent two-agent kinetic

0009-2509/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2013.08.013

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

2

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 1. Molecular reaction structure of P-nitrophenylphosphate and water reacting in presence of acid phosphatase enzyme to form P-nitrophenols and phosphoric acid.

Fig. 2. A simpliﬁed ﬂowchart for the Genetic Algorithm (GA) coupled with Cellular Automata (CA) model system implemented in this study.

interactions of an enzymatic reaction. From the sensitivity analysis study of Dutta et al. (submitted for publication) based on several input parameters, it was realized that there cannot be any speciﬁc mathematical relationship between the input parameters and the model output. It is also not feasible to deﬁne a strict mathematical relation between the chemical kinetic constants and the CA model parameters as the model is an abstract representation of a realistic reaction. This makes repeated trial-and-error interactions of random parameter values the only way to capture speciﬁcally the resulting state of a model system. When there are multiple objective functions of a model to be satisﬁed, useful combinations of parameter values become rarer in the inﬁnite pool of solutions (Konak et al., 2006). If the problem is NP-hard, ﬁnding the optimal solutions by exhaustive search may be too costly or practically impossible. Genetic Algorithm (GA) is a popular meta-heuristic tool frequently used for solving optimization problems (Holland, 1975). Acid phosphatases (APs) are a family of enzymes that are widespread in nature, that non-speciﬁcally catalyze the hydrolysis of orthophosphate monoesters to produce inorganic phosphate and can be found in many animal and plant species (Bull et al., 2002). These enzymes are competitively inhibited by inorganic phosphate (Odds and Hierholzer, 1973). Alvarez (1961) was the ﬁrst to characterize the kinetics and mechanism of the hydrolysis of phosphoric acid esters by potato acid phosphatase. He found that the enzymatic activity depended on three ionizable groups. Based on these ﬁndings, a reaction mechanism was proposed (see Fig. 1). After the nucleophilic attack of one of the ionizable groups on the phosphorous atom of p-nitrophenyl phosphate, the enzyme–substrate complex is formed. After releasing p-nitrophenol, the enzyme–phosphate complex gets hydrolyzed, resulting in the formation of inorganic phosphate and recycled acid phosphatase.

Fig. 3. Reaction steps for the pseudo-artiﬁcial competitive enzyme inhibition reaction representation in the CA model. The symbols indicate Substrate (S), Enzyme (E), Product (P), Inhibitor (I) and Water (W).

In this study, an NSGA-II based Multi-Objective Genetic Algorithm (MOGA) approach is implemented for a stepwise improvement of the CA model. Finally, this approach is validated using the enzymatic hydrolysis of p-nitrophenyl phosphate by acid phosphatase with inhibition reaction.

2. Cellular automata model In accordance with our previous studies, the CA model is framed on a two-dimensional grid. Extended von-Neumann neighborhoods are considered for application of local rules. To represent the boundary as virtually inﬁnite, periodic boundary condition are implemented. This facilitates the boundary molecules to interact with the molecules present in the boundary of opposite direction, as on the surface of toroid. Also, such a neighborhood implementation facilitates the application of probability rules in a combinatorial fashion without any exception due to edge restrictions (Dutta et al., 2011). Each of the reaction compound molecules is considered as a single agent in the model. Each cross section of the grid represents a single lattice

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

node. An agent is allowed to occupy only one lattice node. The types of agents with their percentage occupancy are predeﬁned, and the probability rules are declared as input parameters along with the grid size and iterations. The model can support a maximum of two agent interaction (ﬁrst order and second order reactions) with equal numbers of molecules on either side of any reaction. This means, if a reaction complex forms at any time of the reaction, the model is valid for complexes so formed with two individual agents. The algorithm for rule application remains similar to Dutta et al. (submitted for publication). As earlier, the model is structured on several stochastic rules. There are four basic types of stochastic rules: the probability of Join (PJ) to form complexes, the probability of break (PB) to break a formed complex, the probability of transition (PT) to deﬁne the reaction transition and the free movement probability (PM) which captures the Brownian mixing phenomenon in the model. For a particular agent, each of its neighboring directions is independently examined for feasible events, each of whose probability is obtained in a multiplicative manner. A re-normalization of the directional event probability value is performed. A random direction is then marked for occurrence of the designated event. This makes the model a pure instantaneous model, rather than a memory based model. Each agent is updated once per iteration. Rules are applied in a continuous manner based on the fundamental hypothesis, no two events occur simultaneously.

3. Genetic algorithm for multi-objective optimization Genetic Algorithm is inspired by Darwin's law of survival of ﬁttest (Goldberg, 1989). It mimics the process of natural evolution. It uses operators inspired by biological processes of inheritance, mutation, selection and crossover. The solutions are represented as chromosomes. Each chromosome is composed of values that can either be binary, real or ﬂoating point in representation. Using a set of Objective Functions (OFs), the chromosomes are mapped from variable (genotypic) space to objective (phenotypic) space. A set of solutions deﬁnes the population of a generation. The population evolves as more ﬁt solutions (chromosomes) replaces less ﬁt solutions. This way GA gradually nears to optimal solutions. A multi-objective optimization problem can be stated as a minimization problem as in the following equation: minimize ðf 1 ðxÞ; f 2 ðxÞ; …; f M ðxÞÞ A ℜM where x ¼ ðx1 ; x2 ; …; xn ÞT A ℜn ; and xLi r xi rxUi

ð1Þ

here n is the number of variables. Let u ¼ ðu1 ; u2 ; …; un ÞT A ℜn and v ¼ ðv1 ; v2 ; …; vn ÞT A ℜn be two solutions. Then u dominates v (oru ≻ v), if f i ðuÞ r f i ðvÞ; 8 i and ( j such that f j ðuÞ o f j ðvÞ. Again u and v are said to be nondominated to each other (or ujjv) if neither u≻v nor v≻u. A solution that is not dominated by any other possible

30

100

25

objective function 2

80

objective function 2

3

60

40

20

20 15 10 5

0 0.0000 0.0006 0.0012 0.0018 0.0024 0.0030 0.0036 0.0042 0.0048

0 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014

objective function 1

objective function 1

80 400

70

objective function 2

objective function 2

60 50 40 30 20

300

200

100

10 0 0.0000 0.0004 0.0008 0.0012 0.0016 0.0020 0.0024 0.0028

objective function 1

0

0.00

0.05

0.10

0.15

0.20

0.25

0.30

0.35

objective function 1

Fig. 4. (a) Pareto front at 250th generation using 110 110 grid size obtained from the two objective functions (Eqs. 2a and (2b)), for each member of the population. (Target S conversion ¼23%). (b) Pareto front at 250th generation using 110 110 grid size obtained from the two objective functions (Eqs. (2a) and (2b)), for each member of the population. (Target S conversion ¼18%). (c) Pareto front at 250th generation using 110 110 grid size obtained from the two objective functions (Eqs. (2a) and (2b)). for each member of the population. (Target S conversion ¼10%). (d) Pareto front at 250th generation using 110 110 grid size obtained from the two objective functions (Eqs. (2a) and (2b)). for each member of the population. (Target S conversion¼ 90%).

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

solution is called the Pareto optimal solution. The set of all Pareto optimal solutions is called the Pareto Set (PS). The set of all corresponding objective vectors to PS is called the Pareto front (PF). Given a set of solutions A, a subset of solutions A1 A is obtained such that A1 is a nondominated set of solutions and ðAA1 Þ is a dominated set of solutions. Then A1 becomes the ﬁrst front and nondomination rank of all solutions of A1 is 1. IfðAA1 Þ af, another set of solutions A2 ðAA1 Þ is obtained such that A2 is a nondominated set of solutions. Then A2 becomes the second front and all the solutions of A2 have nondomination rank 2. Similarly an ith front can be deﬁned and therefore, all the solution of ith will have nondomination rank i. Any MOGA tries to ﬁnd a set of solutions which is close to PF, well spread and covers the whole spectrum of PF. Further details on multi objective optimization using genetic algorithms can be found in Deb et al. (2001) and Coello et al. (2007). Several multi-objective genetic algorithms can be found in literature (Zitzler et al., 2001; Deb et al., 2002; Nag and Pal, 2012). For further reviews on MOGA, the interested reader is referred to Coello et al. (1999, 2006) and Konak et al. (2006). NSGA-II (Deb et al., 2002) is an elitist, population based genetic algorithm. The diversity of the population for a particular generation is measured by its crowding distance. Crowding distance can be assigned

to each of the solutions of any nondominated set of solutions. The extreme (boundary) solutions are assigned inﬁnity as crowding distance and for other solutions half the perimeter of its bounding cuboid is set as their crowding distance. It can be conceptualized that in each corner of the cuboid there exists the nearest neighbor according to the corresponding objective axis. The crowding distance assignment procedure (Deb et al., 2002) is shown below: assignCD(SolutionSet P) for each i, P[i].CD ¼0; end for for each objective m, P ¼sortAscending(P, m); P[1].CD ¼INFINITY; P[|P|].CD ¼INFINITY; for i ¼2 to (|P|-1), P[i].CD¼P[i].CD þ(P[|i+1|].m-P[i-1].m) /(P[|P|].m-P[1].m); end for end for end assignCD

350 300

objective function 2

objective function 2

800

600

400

250 200 150 100

200

50 0

0 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 0.0016

0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008

objective function 1

objective function 1

3500 600 3000

objective function 2

objective function 2

500 400 300 200 100

2500 2000 1500 1000 500

0 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0.0040

objective function 1

0 0.00

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

objective function 1

Fig. 5. (a) Pareto front at 250th generation using 310 310 grid size obtained from the two objective functions (Eqs. (2a) and (2b)) for each member of the population. (Target S conversion ¼23%). (b) Pareto front at 250th generation using 310 310 grid size obtained from the two objective functions (Eqs. (2a) and (2b)) for each member of the population. (Target S conversion ¼ 18%). (c) Pareto front at 250th generation using 310 310 grid size obtained from the two objective functions (Eqs. (2a) and (2b)) for each member of the population. (Target S conversion ¼10%). (d) Pareto front at 250th generation using 310 310 grid size obtained from the two objective functions (Eqs. (2a) and (2b)) for each member of the population. (Target S conversion ¼90%).

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4. NSGA-II algorithm Step-I (Population initialization): randomly initializes N number of solutions, where N is the population size. Step-II (Generate N number of offspring): a mating pool is created using a binary tournament selection where a two-tier

Fig. 6. (a) Molecular structure of inhibitor sodium phosphate monohydrate. (b) Intermediate steps of the enzymatic reaction (S: P-nitrophenylphosphate, E: acid phosphatase enzyme, A: Pi-enzyme complex, B: P-nitrophenol, W: water, P: Phosphate ion and I: inhibitor) used for model validation in the study.

5

ﬁtness mechanism is used. The solution with lower rank is preferred. If the solutions have the same nondomination rank, then the solution with higher crowding distance is the winner. Then, Simulated Binary Crossover (SBX) (Deb and Agrawal, 1995) is performed on the mating pool to create N offspring. Polynomial mutation (Deb and Tiwari, 2006) is followed by crossover to mutate the offspring. Next, all the mutated offspring are evaluated. Thus a set of N offspring is generated. Step-III (Modify the population): the offspring and the parent population are added together to create a union of size 2N. After that, from these 2N solutions, the best N solutions are selected to form the next generation. In this selection procedure, NSGA-II uses a two-tier ﬁtness mechanism. Lower nondomination rank is given the higher priority and kept in the ﬁrst tier. Crowding distance is used in the second tier to differentiate the quality of the solutions having same rank. To frontify, NSGA-II uses fast-non-dominated-sort (Deb et al., 2002) procedure. Step-IV (Terminate the algorithm): if the termination condition (may be a maximum number of function evaluations or generations) is satisﬁed, perform a fast-non-dominated-sort and

Table 1a Optimized results of the coupled CA–GA simulation using initial concentrations (Substrate 1 (S)¼ 43.12%, Substrate 2 (W) ¼1.31%, Inhibitor (I)¼ 21.56%, Enzyme (E) ¼ 0.005% and Void (_) ¼ 34%) for a 150 150 CA grid size, 5000 CA iterations, 100 GA population size and 250 GA generations of simulation. Final target Product (P) ¼1.30%. The chemical reaction time is 30 min. Simulation time with 36 parallel threads took approximately 6 h. P_Join_ES

P_Join_AB

P_Join_AW

P_Join_EI

P_Break_ES

P_Break_AB

P_Break_EP

P_Break_EI

P_trans_ESAB

P_trans_ABES

P_Trans_AWEP

Iterations

%P

0.702837 0.653945 0.613252

0.121707 0.177088 0.370823

0.781171 0.817983 0.790154

0.029324 0.887521 0.797484

0.133784 0.393030 0.349440

0.590446 0.868474 0.600341

0.622463 0.839000 0.358055

0.079896 0.922149 0.953598

0.978648 0.636805 0.646044

0.312356 0.131928 0.269840

0.990690 0.990050 0.991711

163,552 183,957 215,268

1.30 1.30 1.30

Table 1b Optimized results of the coupled CA–GA simulation using initial concentrations (Substrate 1 (S) ¼38.39%, Substrate 2 (W) ¼8.41%, Inhibitor (I) ¼19.19%, Enzyme (E) ¼ 0.002823% and Void (_) ¼34%) for a 150 150 CA grid size, 5000 CA iterations, 100 GA population size and 250 GA generations of simulation. Final target Product (P) ¼ 8.41%. The chemical reaction time is 2.5 min. Simulation time with 36 parallel threads took approximately 6 h. P_Join_ES

P_Join_AB

P_Join_AW

P_Join_EI

P_Break_ES

P_Break_AB

P_Break_EP

P_Break_EI

P_trans_ESAB

P_trans_ABES

P_Trans_AWEP

Iterations

%P

0.031686 0.117018 0.985828 0.034687 0.956729

0.044148 0.037198 0.681940 0.442827 0.005193

0.890304 0.984235 0.979454 0.979739 0.986940

0.117166 0.160812 0.492411 0.163103 0.457661

0.524927 0.520963 0.649391 0.557778 0.651398

0.991674 0.868183 0.719708 0.865898 0.945841

0.711479 0.311518 0.619185 0.311518 0.619185

0.710109 0.866367 0.740343 0.866367 0.740930

0.531081 0.959528 0.950976 0.948876 0.968777

0.001208 0.028873 0.012735 0.013018 0.024097

0.994324 0.995971 0.979529 0.995933 0.979792

272,196 400,000 289,698 321,388 400000

8.41 8.39 8.41 8.40 8.39

Table 1c Optimized results of the coupled CA–GA simulation using initial concentrations (Substrate 1 (S) ¼ 35.31%, Substrate 2 (W) ¼13.04%, Inhibitor (I) ¼17.64%, Enzyme (E) ¼ 0.002596% and Void (_) ¼34%) for a 150 150 CA grid size, 5000 CA iterations, 100 GA population size and 250 GA generations of simulation. Final target Product (P) ¼ 13.04%. The chemical reaction time is 5 min. Simulation time with 36 parallel threads took approximately 6 h. P_Join_ES

P_Join_AB

P_Join_AW P_Join_EI

P_Break_ES P_Break_AB P_Break_EP P_Break_EI P_trans_ESAB P_trans_ABES P_Trans_AWEP

Iterations % P

0.477159 0.465103 0.774559 0.765923 0.767359

0.034581 0.045059 0.436236 0.350941 0.033105

0.773444 0.805223 0.805223 0.834075 0.688106

0.086241 0.078904 0.078919 0.531113 0.153293

400,000 400,000 400,000 400,000 400,000

0.029921 0.268620 0.268620 0.633454 0.274920

0.884724 0.905332 0.989605 0.989776 0.979561

0.936365 0.890624 0.958038 0.991090 0.981325

0.943321 0.888123 0.888123 0.810838 0.329886

0.922237 0.857370 0.857370 0.947486 0.945152

0.008024 0.020236 0.020236 0.017457 0.017488

0.989975 0.988668 0.988668 0.960751 0.986868

11.08 13.00 13.01 12.98 13.03

Table 1d Optimized results of the coupled CA–GA simulation using initial concentrations (Substrate 1 (S) ¼27.81%, Substrate 2 (W) ¼24.27%, Inhibitor (I) ¼13.90%, Enzyme (E) ¼ 0.0204% and Void (_) ¼34%) for a 150 150 CA grid size, 5000 CA iterations, 100 GA population size and 250 GA generations of simulation. Final target Product (P) ¼24.27%. The chemical reaction time is 2.5 min. Simulation time with 36 parallel threads took approximately 6 h. P_Join_ES

P_Join_AB

P_Join_AW P_Join_EI

0.816042 0.993704

0.031226 0.103822

0.811154 0.811589

P_Break_ES P_Break_AB P_Break_EP P_Break_EI P_trans_ESAB P_trans_ABES P_Trans_AWEP

0.258073 0.006672 0.846079 0.007089

0.961697 0.999228

0.921515 0.948470

0.998066 0.915764

0.999853 0.999829

0.000574 0.002215

0.994197 0.995784

Iterations % P 400,000 400,000

24.26 24.25

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

return the nondominated solutions of the archive, else go to Step-II.

5. Model formulation The CA model implemented in this study acts as a black box for NSGAII. The input parameters are selected by the GA, which is then used to simulate the CA model. The simulation output is then used to calculate the ﬁtness functions. Finally, new sets of input parameter values are obtained after implementing the genetic operators of the GA. The model iterates until the ﬁnal targeted generation count is reached. The following parameters are considered as target optimization parameters: the joining probabilities between enzyme and substrate (Join (ES)), between enzyme and inhibitor (Join (EI)), between enzyme and product (Join (EP)), between water and substrate (Join (WS)) and between two water agents (Join (WW)), the breaking probabilities between enzyme and substrate (Break (ES)), between enzyme and inhibitor (Break (EI)), between enzyme and product (Break (EP)), between substrate and water (Break (WS)) and between two water agents (Break (WW)), the transition probabilities from enzyme–substrate to enzyme– product (Trans (ESEP)) and from enzyme–product to enzyme– substrate (Trans (EPES)). Due to the probabilistic nature of these parameters, the values are varied within the range of 1.0E-10–9.99E-10 (correct to ten decimal places). Each of these parameters can have any real value within the range.

The optimization approach used in this study is based on mainly three objective functions (OFs). The ﬁrst OF (Eq. (2a)) is the summation of all the probabilistic parameter values. The second OF (Eq. (2b)) is the absolute difference of the product concentration between the target and the simulation output, thereby to obtain the model output product concentration as close as possible to the target product concentration. A third OF (Eq. (2c)) has been tested during model validation work and is used to replace the ﬁrst OF. It is deﬁned to measure the initial speed of reaction for product formation. This initial speed of the reaction is deﬁned as ﬁrst few n iterations of the model. The iteration number can be deﬁned as the initial simulation settings. The ﬁrst two functions are targeted to be minimized using the NSGA-II algorithm. For the third OF, whose target is to maximize in the positive axis, the effect is achieved by considering the negative of the value. The minimization of the summation of the parameter values (ﬁrst OF) facilitates the optimizer to optimize the function while keeping the individual parameter values to minimum. Since the initial state of the CA model considers product concentration to be zero, these probabilistic parameter values have to increase from zero towards unity to reach the target product concentration. As such the summation value will always have an increasing tendency so that the model output product concentration can reach nearer to the target product concentration. This makes the ﬁrst OF to have a mutually inverse relationship with the expected property. Furthermore, the second OF facilitates increase of ﬁnal product formation and the third OF targets to increase the initial speed of the reaction so that more products are formed with fewer iterations. The ﬂowchart in Fig. 2 shows the logical ﬂow structure

Fig. 7. (a–d): Substrate concentration proﬁles obtained from the simulation runs using optimized parameter sets as mentioned in Tables (1a–1d) respectively. The x-axis represents iteration number and the y-axis represents the substrate count at its corresponding iteration step. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

of the GA–CA coupled approach. ! objective f 1 ¼ min

12

ð2aÞ

∑ Pi

i¼1

objective f 2 ¼ minjP target P model j

ð2bÞ

objective f 3 ¼ maxðP target P n Þ

ð2cÞ

The CA program is developed in the C language. The Java coded open source jMETAL Optimization package (Durillo and Nebro, 2011) is used for the NSGA II algorithm. The communication between the GA and the CA occurs through ﬁle based I/O operation. This helps the simulation design to be robust as the simulations can be restarted from the ﬁnal terminated GA generation by using the status storage ﬁles.

6. Results and discussions Initially, some pseudo real reaction sets are considered in order to study the capability of the model. A competitive enzyme inhibition reaction scheme is represented for all the simulation runs as shown in Fig. 3. Such a reaction structure is modeled so that one can study the ability of the GA optimization. The representation considers substrate–water afﬁnity and afﬁnity of two water molecules among each other. This also illustrates the possibility of parallel reactions. The reaction scheme, in total, has twelve probability rule parameters. The optimization study is performed in two grid sizes namely, 110 110 and 310 310 so as to compare the results. The initial parameter settings are set for four different target Substrate Conversion to Product (ScP): 10%, 18%, 23% and 90%. The initial

7

concentration of water, substrate, enzyme, inhibitor and product are ﬁxed to 60%, 4%, 1%, 1% and 0% respectively. 34% of void is maintained as per the molecular theory of intermediate space. The total number of iteration steps for each CA is ﬁxed to 3000. Each optimization simulation is run for 250 generations, with 100 as the population size of each generation. Each of the single threaded optimization simulations for a 310 310 grid took approximately 22 days, while each of the optimization simulations for a 110 110 grid took approximately 4 days for completion, with the HPC Flemish computing facility. Fig. 4(a–d) shows the characteristic Pareto front for the 110 110 grid for the corresponding ScP. The Pareto fronts for the 310 310 grid are shown in Fig. 5(a–d). There can be many combinations of probability parameter values that can be considered for simulation to reach a single ScP within considerable SD. Hence each optimization simulation results in multiple sets of near-optimal parameter values. But all these result sets are not realistic as some may have highly optimized ﬁrst OF which may lead to poor second OF. To ﬁlter out such anomalies, the Two sample Z-Test method (see Eq. (3)) is implemented on the simulation results sets, after calculating mean and standard deviation over 100 runs using the following relation: jM E M i j Z ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ s2E s2i

ð3Þ

where ME is expected mean P count, Mi is obtained mean P count, sE and si are the standard deviations of ME and Mi respectively. A parameter set is considered acceptable only if Z value is less than 2. The acceptable detailed result sets and the corresponding probability parameter values are listed in Appendix (see Tables A1–A8).

Table A1 Grid size: 110 110. Target % conversion of S to P: 10%. The mean is calculated over 100 runs for each of the parameter sets. Standard deviation is calculated over the mean. % S to P shows the percentage of mean S converted to P after completion of simulation. The probability values are shown correct up to 12 decimal places. This optimization run gave 6 optimized parameter sets. Mean [P]

Standard deviation

% S to P

21.97

4.250858499

4.539256198

22.26

4.384591938

4.599173554

28.49

5.238734908

5.886363636

30.06

4.956008495

6.210743802

31.25

4.660981385

6.45661157

35.83

5.915319949

7.402892562

Mean [P]

Standard deviation

% S to P

21.97

4.250858499

4.539256198

22.26

4.384591938

4.599173554

28.49

5.238734908

5.886363636

30.06

4.956008495

6.210743802

31.25

4.660981385

6.45661157

35.83

5.915319949

7.402892562

Probability values P_Join_ES

P_Join_EI

P_Join_EP

P_Join_WW

P_Join_WS

P_Break_ES

1.000602917776 E 04 1.000707755436 E 04 1.000651238129 E 04 1.000485570402 E 04 1.000485570402 E 04 1.011290550471 E 04

1.007932744715 E 04 1.008477365938 E 04 1.008563500405 E 04 1.008017845003 E 04 1.008017845003 E 04 1.008399050059 E 04

1.001266200158 E 04 1.008864737424 E 04 1.008959081391 E 04 1.008963274401 E 04 1.008963274401 E 04 1.003920087651 E 04

1.001520962471 E 04 1.000807130551 E 04 1.001500379831 E 04 1.000797799512 E 04 1.000797799512 E 04 1.000333779073 E 04

1.091070808353 E 04 1.000182733831 E 04 1.000182238827 E 04 1.000020350681 E 04 1.000182238827 E 04 1.000879940322 E 04

1.001391492187 E 04 1.000014507478 E 04 1.000888000415 E 04 1.000888958020 E 04 1.000888958020 E 04 1.001053294020 E 04

P_Break_EI

P_Break_EP

P_Break_WW

P_BreakWS

P_Trans_ESEP

P_Trans_EPES

1.036741468500 E 04 1.000222682473 E 04 1.001872655570 E 04 1.000089641649 E 04 1.000089641649 E 04 1.027779078866 E 04

1.046394902820 E 04 1.001406194346 E 04 1.011912091564 E 04 1.000197016245 E 04 1.000197016245 E 04 1.005681258140 E 04

1.004419952226 E 04 1.004742423281 E 04 1.000373243938 E 04 1.001214997529 E 04 1.001214997529 E 04 1.004490776696 E 04

1.002645259265 E 04 1.000286603630 E 04 1.002790531377 E 04 1.010796705201 E 04 1.010796705201 E 04 1.002284890614 E 04

2.066305374693 E 04 1.889213316691 E 04 2.646054958676 E 04 2.882219143478 E 04 2.882219143478 E 04 3.717183520624 E 04

1.007793267166 E 04 1.001391520289 E 04 1.000192533237 E 04 1.001074716142 E 04 1.001788998350 E 04 1.001309340156 E 04

Probability values

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

7. Model validation To represent any kinetic reaction through a CA model, some informations are needed a priori such as the intermediate steps of

the reaction, the initial concentration of the reacting agents and the ﬁnal concentration of the target product. Since the CA model considers one node to one agent representation, the concentration of the reaction agents can be converted to % occupancy based on

Table A2 Grid size: 310 310. Target % conversion of S to P: 10%. The mean is calculated over 100 runs for each of the parameter sets. Standard deviation is calculated over the mean. % S to P shows the percentage of mean S converted to P after completion of simulation. The probability values are shown correct up to 12 decimal places. This optimization run gave 15 optimized parameter sets. Mean [P]

Standard deviation

% S to P

319.38

14.62098948

8.308532778

319.73

15.77035641

8.317637877

320.74

14.81578465

8.343912591

320.83

16.71451012

8.346253902

323.2

17.61484681

8.407908429

324.56

14.94793658

8.443288241

328.04

16.20320457

8.533818939

328.05

19.16561262

8.534079084

333.14

15.26865145

8.666493236

334.52

16.31351921

8.70239334

335.34

18.21932596

8.723725286

357.2

16.82230362

9.292403746

358.17

17.14528236

9.317637877

359.1

17.21433511

9.341831426

363.61

19.23690424

9.459157128

Mean [P]

Standard deviation

% S to P

319.38

14.62098948

8.308532778

319.73

15.77035641

8.317637877

320.74

14.81578465

8.343912591

320.83

16.71451012

8.346253902

323.2

17.61484681

8.407908429

324.56

14.94793658

8.443288241

328.04

16.20320457

8.533818939

328.05

19.16561262

8.534079084

333.14

15.26865145

8.666493236

334.52

16.31351921

8.70239334

335.34

18.21932596

8.723725286

357.2

16.82230362

9.292403746

358.17

17.14528236

9.317637877

359.1

17.21433511

9.341831426

363.61

19.23690424

9.459157128

Probability values P_Join_ES

P_Join_EI

P_Join_EP

P_Join_WW

P_Join_WS

P_Break_ES

1.010116078510 E 04 1.010108491114 E 04 1.010100844735 E 04 1.001324819845 E 04 1.009674694182 E 04 1.009373305854 E 04 1.008335277886 E 04 1.009418953694 E 04 1.010116078510 E 04 1.010116078510 E 04 1.010116078510 E 04 1.010023360548 E 04 1.009376335052 E 04 1.009376335052 E 04 1.008921580439 E 04

1.015615254564 E 04 1.015606721121 E 04 1.003130306908 E 04 1.001375804160 E 04 1.015606721121 E 04 1.003621417207 E 04 1.002114602378 E 04 1.003621417207 E 04 1.015219480414 E 04 1.015219480414 E 04 1.015219480414 E 04 1.000926981735 E 04 1.015690501595 E 04 1.015690501595 E 04 1.002572808875 E 04

1.018313923547 E 04 1.018202818464 E 04 1.012885464358 E 04 1.012885464358 E 04 1.018310155342 E 04 1.010604910860 E 04 1.015290033590 E 04 1.010644848077 E 04 1.013594603161 E 04 1.057941770484 E 04 1.013594603161 E 04 1.033119394146 E 04 1.018344282947 E 04 1.017843941959 E 04 1.005468108625 E 04

1.001603873841 E 04 1.001808762122 E 04 1.000888230856 E 04 1.000602531220 E 04 1.001808762122 E 04 1.001851993363 E 04 1.000785250360 E 04 1.001836759821 E 04 1.008269541408 E 04 1.008289686581 E 04 1.008289686581 E 04 1.001120798044 E 04 1.050165823899 E 04 1.050165823899 E 04 1.022771546602 E 04

1.009314295535 E 04 1.004729042804 E 04 1.000133492519 E 04 1.000133492519 E 04 1.009322369088 E 04 1.005494086382 E 04 1.001087815211 E 04 1.005494086382 E 04 1.004781002502 E 04 1.004781002502 E 04 1.004781002502 E 04 1.008794096407 E 04 1.009322369088 E 04 1.002744695534 E 04 1.004499064361 E 04

1.007026244194 E 04 1.007026927431 E 04 1.007290993068 E 04 1.007291474909 E 04 1.007026780311 E 04 1.006726876782 E 04 1.006854290852 E 04 1.006726876782 E 04 1.007034848414 E 04 1.007034848414 E 04 1.007034848414 E 04 1.011803363491 E 04 1.006758023702 E 04 1.006758023702 E 04 1.012523544310 E 04

Probability values P_Break_EI

P_Break_EP

P_Break_WW

P_BreakWS

P_Trans_ESEP

P_Trans_EPES

1.019035715092 E 04 1.019332455562 E 04 1.036055127644 E 04 1.021926197124 E 04 1.020161047083 E 04 1.035394766180 E 04 1.038342655211 E 04 1.000734808738 E 04 1.019188632769 E 04 1.019311979509 E 04 1.019035715092 E 04 1.003251853000 E 04 1.007613238847 E 04 1.007613238847 E 04 1.035290524325 E 04

1.021042514323 E 04 1.021315405996 E 04 1.021661249196 E 04 1.004446447096 E 04 1.033810445777 E 04 1.017109318253 E 04 1.021361317088 E 04 1.018549289414 E 04 1.021045348193 E 04 1.021045348193 E 04 1.021045348193 E 04 1.033756835919 E 04 1.021452160614 E 04 1.021452160614 E 04 1.015393643861 E 04

1.023445719230 E 04 1.023442640595 E 04 1.006571340831 E 04 1.006571340831 E 04 1.023443868272 E 04 1.022961444836 E 04 1.036084692034 E 04 1.004721492186 E 04 1.017039762712 E 04 1.023332054088 E 04 1.023332054088 E 04 1.030809078709 E 04 1.030981687059 E 04 1.030981687059 E 04 1.027995994903 E 04

1.005045526915 E 04 1.007549524855 E 04 1.000722908252 E 04 1.000049620473 E 04 1.003234766233 E 04 1.000031218509 E 04 1.019427047391 E 04 1.000031218509 E 04 1.005045526915 E 04 1.007734501307 E 04 1.005045526915 E 04 1.000149708747 E 04 1.002804331523 E 04 1.002804331523 E 04 1.014181178972 E 04

4.432767952952 E 04 4.442083292880 E 04 4.488917161322 E 04 4.508197059966 E 04 4.442083292880 E 04 4.657091409030 E 04 4.660700508197 E 04 4.653476215159 E 04 4.772713665644 E 04 4.852064830617 E 04 4.852064830617 E 04 5.523795109237 E 04 5.429341782261 E 04 5.478424800058 E 04 5.610445209127 E 04

1.003461363248 E 04 1.003460922417 E 04 1.002310736193 E 04 1.002310736193 E 04 1.003460922417 E 04 1.006308254705 E 04 1.003913329033 E 04 1.006308254705 E 04 1.007231072831 E 04 1.003724099618 E 04 1.003724099618 E 04 1.002496418067 E 04 1.003645351088 E 04 1.003724718207 E 04 1.003639501466 E 04

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

the calculation of Avogadro number. The void consideration of 34% on the CA grid is kept constant in all representations. Speciﬁcally for enzymatic reaction data, it is observed that the concentration scales of the reacting substrate and that of the participating enzyme are widely different. Such difference arises from the reusability and high activity properties of the enzyme. To be able to represent such a broad scale of concentrations, the size of the CA grid requires to be increased proportionally. But a larger CA grid will require signiﬁcantly more computational time. Hence, to keep the grid within manageable limit, the minimal concentration representation format is followed. The lowest-concentration agent will have the least possible count in the CA grid, possibly just one. But such a representation can cause reaction completion time reach inﬁnity. After an initial exponential phase of the reaction the number of forward reacting agents reduce. This causes the overall probability of feasible neighborhood interaction events among the agents become unpredictably high. The randomness of free movement increases thereby delaying the expected forward reaction. This can be observed from the fact that the simulation takes only 83,611 iterations for the ﬁrst 98.64% of substrate conversion, whereas it takes as much as 131,657 iterations for the next 1% of substrate conversion. Hence during optimization, the number of iterations required for automata evolution can be set. To choose the number of iterations for a particular reaction structure, sample runs of the CA model are performed with sufﬁciently high number of iterations until the reaction equilibrium is reached, while the

9

values of the probability rules are changed in these simulations. Further while analyzing the concentration curves, additional iterations are considered so that it is well beyond the completion or attainment of reaction equilibrium. In this way, the resulting optimized parameter value sets can be used to simulate the CA model that can reach the target ScP provided. However, it must be noted that a large number of CA iterations may be required to reach the ﬁnal targeted ScP due to the reduction in reaction speed once reaction completion/equilibrium is attained. As mentioned earlier, the enzymatic hydrolysis of pnitrophenyl phosphate by acid phosphatase with inhibition reaction is chosen for model validation. In the experiment, acid phosphatase (AP) from potato (P3752, Sigma Aldrich) is used. The enzymatic function in the production, transport and recycling of phosphate is critical for the metabolic and energy transduction processes of a living cell. In the ﬁrst reaction, the enzyme– substrate complex is formed by a nucleophilic attack on the phosphorous atom of the substrate. Subsequently, p-nitrophenol is released. Next, the enzyme–Pi complex is hydrolyzed to yield phosphate ions and recycled enzyme. Sodium phosphate monohydrate (NaH2PO4 H2O) (Fig. 6a) is used as an inhibitor in the reaction. The intermediate steps of the reaction as represented by the CA–GA model are shown in Fig. 6b. Four sets of experimental data with varying reactant and enzyme concentrations for different reaction times are collected for the validation. The concentration of the enzyme is proportionally less by a large factor of ten for

Table A3 Grid size: 110 110. Target % conversion of S to P: 18%. The mean is calculated over 100 runs for each of the parameter sets. Standard deviation is calculated over the mean. % S to P shows the percentage of mean S converted to P after completion of simulation. The probability values are shown correct up to 12 decimal places. This optimization run gave 8 optimized parameter sets. Mean [P]

Standard deviation

% S to P

50.2

6.237715199

10.37190083

50.23

6.524011594

10.37809917

51.65

7.061597667

10.6714876

58.62

6.869313994

12.11157025

58.81

6.064793251

12.15082645

60.23

7.832662769

12.44421488

63.43

7.210058911

13.1053719

69.28

7.562640766

14.31404959

Mean [P]

Standard deviation

% S to P

50.2

6.237715199

10.37190083

50.23

6.524011594

10.37809917

51.65

7.061597667

10.6714876

58.62

6.869313994

12.11157025

58.81

6.064793251

12.15082645

60.23

7.832662769

12.44421488

63.43

7.210058911

13.1053719

69.28

7.562640766

14.31404959

Probability values P_Join_ES

P_Join_EI

P_Join_EP

P_Join_WW

P_Join_WS

P_Break_ES

1.048094712545 E 04 1.001244347383 E 04 1.053479852615 E 04 1.001872665734 E 04 1.053479852615 E 04 1.143517759875 E 04 2.254269386759 E 04 1.054961553482 E 04

1.001505475557 E 04 1.001559755990 E 04 1.001046551604 E 04 1.001597128727 E 04 1.001047752925 E 04 1.001546537871 E 04 1.249581571941 E 04 1.001046551604 E 04

1.008920356947 E 04 1.001224952003 E 04 1.003924042678 E 04 1.007029854363 E 04 1.003826491327 E 04 1.130420962239 E 04 1.123985437302 E 04 1.004794125149 E 04

1.021560028520 E 04 1.004675642990 E 04 1.006872681835 E 04 1.023062702913 E 04 1.006872681835 E 04 1.114382387083 E 04 1.135813364430 E 04 1.006872681835 E 04

1.026789654284 E 04 1.111230565156 E 04 1.043339722100 E 04 1.062597984631 E 04 1.043339722100 E 04 1.119155661982 E 04 1.066708078483 E 04 1.043339722100 E 04

1.014229199687 E 04 1.000551083323 E 04 1.005184187319 E 04 1.009222413425 E 04 1.008987870839 E 04 1.001695234555 E 04 1.006245812753 E 04 1.000668824351 E 04

P_Break_EI

P_Break_EP

P_Break_WW

P_BreakWS

P_Trans_ESEP

P_Trans_EPES

1.099928433096 E 04 1.028091644992 E 04 1.006469305580 E 04 1.011958030237 E 04 1.006469305580 E 04 1.263394826444 E 04 1.653470872448 E 04 1.004285804416 E 04

1.004381784905 E 04 1.000442974830 E 04 1.002861241813 E 04 1.001081163835 E 04 1.003069106566 E 04 1.119805201812 E 04 1.003584587830 E 04 1.002465316386 E 04

1.156767553391 E 04 1.000229984753 E 04 1.137723153062 E 04 1.000520454660 E 04 1.137170899909 E 04 1.089570342539 E 04 1.145986392025 E 04 1.138187025268 E 04

1.010398023701 E 04 1.011931635005 E 04 1.145341489468 E 04 1.014406544194 E 04 1.150571751737 E 04 1.014372402018 E 04 1.023501909793 E 04 1.010145680744 E 04

7.016775863385 E 04 7.023906026078 E 04 7.207797739175 E 04 9.428807124774 E 04 9.711369160627 E 04 1.064254441294 E-03 1.266878814938 E-03 1.659525058204 E-03

1.132892327710 E 04 1.129454482250 E 04 1.028957045455 E 04 1.021986301730 E 04 1.029021208279 E 04 1.140716701943 E 04 1.220892253579 E 04 1.028957045455 E 04

Probability values

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

which the minimal agent representation format is followed, as described earlier. The CA iteration is ﬁxed to 5000. Each of the GA optimization runs for 250 generations with population size 100. Eqs. (2b) and (2c) are the two optimization functions used for the validation scheme. The concentration gradient for the initial speed of reaction (for Eq. (2c)) is measured for the ﬁrst 10,000 CA iterations. The results of the GA optimization are shown in Tables 1a–1d. Each element of the optimization set has a varying number of resultant optimized parameter sets. Due to the stochastic nature, the number of obtained optimized parameter sets may vary for different runs and for different targets of optimization. The obtained optimal parameters are then used to run independent simulations, each for a considerable number of CA iterations steps. The concentration graphs are plotted along the number of iterations for corresponding reaction set, as shown in Fig. 7(a–d). The general characteristics of these graphs match well with the concentration proﬁles. Although most of the curves are overlapping, at times the characteristics can be drastically affected due to stochastic noise, as seen in Fig. 7c (red curve, line1). Such a noise can be minimized only if the observations are made over a mean of several runs. However due to extremely high computational cost, this approach was avoided. The product yield (%P) obtained from the simulation tables and the simulation curves indicates the validity of the CA–GA optimization tool. The CPU time for each simulation is approximately 6 h. The simulations are run using a high performance compute cloud available on credit basis from Amazon Web Service (AWS). Each simulation is enabled

to run with 36 parallel threads in dedicated processing cores available through AWS hardware level virtualization.

8. Conclusion Cellular automata can represent a discrete system with a simple yet abstract representation of a natural process governed by probability rules. However, an abstract representation makes it highly inadequate to have a mathematical relationship between the process parameters, so that parameter estimation by trial and error is the only choice left to synchronize the results. In the current study, it is successfully shown that genetic algorithms can be used to optimize the parameter sets of a kinetic model. The validation study using enzymatic hydrolysis of disodium phosphate by acid phosphatase with inhibition reaction aptly describes the usefulness of the CA model when used with a GA optimizer to represent (bio)chemical reaction in silico. A practical implication of this study is that the stochastic CA model can represent each and every intermediate step of the kinetic reaction, which possibly cannot be observed in a deterministic model due to its limitations. Further it can be used to represent complex biological reactions and the simulation results can be analyzed without actually performing in vitro. It is also shown that an extensive need for computational resources for discrete modeling can be met using high performance cloud clusters and farm computational clusters, which are available at a reasonable price.

Table A4 Grid size: 310 310. Target % conversion of S to P: 18%. The mean is calculated over 100 runs for each of the parameter sets. Standard deviation is calculated over the mean. % S to P shows the percentage of mean S converted to P after completion of simulation. The probability values are shown correct up to 12 decimal places. This optimization run gave 8 optimized parameter sets. Mean [P]

Standard deviation

% S to P

Probability values P_Join_ES

613.56

23.52575206

15.96149844

614.74

22.19137151

15.99219563

615.57

22.92731378

16.01378772

615.87

23.0496412

16.02159209

616.24

21.74141883

16.03121748

623.28

20.17052555

16.21436004

632.65

22.92912884

16.45811655

639.22

23.54260269

16.62903226

Mean [P]

Standard deviation

% S to P

613.56

23.52575206

15.96149844

614.74

22.19137151

15.99219563

615.57

22.92731378

16.01378772

615.87

23.0496412

16.02159209

616.24

21.74141883

16.03121748

623.28

20.17052555

16.21436004

632.65

22.92912884

16.45811655

639.22

23.54260269

16.62903226

1.019545281462 E 04 1.208458855382 E 04 1.007602216705 E 04 1.019545281462 E 04 1.006979465550 E 04 1.013234359529 E 04 1.213749989209 E 04 1.208458855382 E 04 Probability values

P_Join_EI

P_Join_EP

P_Join_WW

P_Join_WS

P_Break_ES

1.104395051114 E 04 1.250517100276 E 04 1.136848865553 E 04 1.141757495236 E 04 1.137110834471 E 04 1.252313439663 E 04 1.249206318580 E 04 1.250517100276 E 04

1.018247258238 E 04 1.018623295369 E 04 1.018283582553 E 04 1.018254003529 E 04 1.018254003529 E 04 1.018079236632 E 04 1.112333149000 E 04 1.116361096101 E 04

1.002323794526 E 04 1.000303132528 E 04 1.006868214435 E 04 1.002323794526 E 04 1.006872557415 E 04 1.002323794526 E 04 1.013726574571 E 04 1.000275020031 E 04

1.008198024590 E 04 1.000385219250 E 04 1.005152126541 E 04 1.011304568807 E 04 1.011304568807 E 04 1.553518063714 E 04 1.278605438845 E 04 1.548351567803 E 04

1.086843790134 E 04 1.090639708781 E 04 1.003770737595 E 04 1.086698518660 E 04 1.086698518660 E 04 1.086698518660 E 04 1.106510291229 E 04 1.090639708781 E 04

P_Break_EI

P_Break_EP

P_Break_WW

P_BreakWS

P_Trans_ESEP

P_Trans_EPES

1.008945422625 E 04 1.416309274615 E 04 1.018564936277 E 04 1.008945422625 E 04 1.018564936277 E 04 1.393180927301 E 04 1.007916449738 E 04 1.410891812727 E 04

1.036427123396 E 04 1.192955659244 E 04 1.000001993533 E 04 1.036427123396 E 04 1.000914622733 E 04 1.184038323889 E 04 1.188544422288 E 04 1.192955659244 E 04

3.279509642706 E 04 3.033196981411 E 04 3.210513228082 E 04 3.279509642706 E 04 3.279509642706 E 04 3.287376012639 E 04 4.404104836872 E 04 4.346235894080 E 04

1.014181160988 E 04 1.186516396678 E 04 1.009950712949 E 04 1.018718153486 E 04 1.009950712949 E 04 1.089535892753 E 04 1.023345522340 E 04 1.186516396678 E 04

1.949823836323E-03

1.038727299869 E 04 1.038994559704 E 04 1.038727299869 E 04 1.038727299869 E 04 1.038727299869 E 04 1.019247662717 E 04 1.019322322930 E 04 1.019322322930 E 04

1.955588129375E-03 1.949823836323E-03 1.949823836323E-03 1.949823836323E-03 1.949823836323E-03 1.955588129375E-03 1.955588129375E-03

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Acknowledgment SK and AD would like to acknowledge Advait Apte, Arizona State University (USA) for his valuable inputs and discussions on

11

the cellular automata model. AD would like to thank Paul Vandecruys, Laboratory of Molecular Cell Biology, Institute of Botany and Microbiology, KU Leuven (Belgium) for his help with

Table A5 Grid size: 110 110. Target % conversion of S to P: 23%. The mean is calculated over 100 runs for each of the parameter sets. Standard deviation is calculated over the mean. % S to P shows the percentage of mean S converted to P after completion of simulation. The probability values are shown correct up to 12 decimal places. This optimization run gave 15 optimized parameter sets. Mean [P]

Standard deviation

% S to P

64

6.999278462

13.2231405

64.65

6.507765369

13.35743802

65.49

6.818713111

13.53099174

66.17

7.36227759

13.6714876

66.6

7.51295178

13.76033058

67.51

6.753218462

13.94834711

68.59

7.96215163

14.1714876

68.7

8.022065529

14.19421488

71.96

8.501835809

14.8677686

72.23

7.069431992

14.92355372

75.45

7.658592586

15.58884298

77.08

8.049945103

15.92561983

80.66

8.065664851

16.66528926

81.65

8.035747154

16.86983471

86.83

9.850452484

17.94008264

Mean [P]

Standard deviation

% S to P

64

6.999278462

13.2231405

64.65

6.507765369

13.35743802

65.49

6.818713111

13.53099174

66.17

7.36227759

13.6714876

66.6

7.51295178

13.76033058

67.51

6.753218462

13.94834711

68.59

7.96215163

14.1714876

68.7

8.022065529

14.19421488

71.96

8.501835809

14.8677686

72.23

7.069431992

14.92355372

75.45

7.658592586

15.58884298

77.08

8.049945103

15.92561983

80.66

8.065664851

16.66528926

81.65

8.035747154

16.86983471

86.83

9.850452484

17.94008264

Probability values P_Join_ES

P_Join_EI

P_Join_EP

P_Join_WW

P_Join_WS

P_Break_ES

1.051787775023 E 04 1.018315667165 E 04 1.022879761275 E 04 1.055845332059 E 04 1.018315667165 E 04 1.022991309271 E 04 1.013672387770 E 04 1.276648168719 E 04 1.276626557766 E 04 1.279310133556 E 04 1.275724029806 E 04 1.276648168719 E 04 1.016965800795 E 04 1.126431018368 E 04 1.840930741906 E 04

1.002749689329 E 04 1.011717193751 E 04 1.005678948279 E 04 1.039815327125 E 04 1.011717193751 E 04 1.056288642630 E 04 1.225827910553 E 04 1.040586827613 E 04 1.040586827613 E 04 1.043382399437 E 04 1.204096701189 E 04 1.040586827613 E 04 1.268423478397 E 04 1.197410555395 E 04 1.226871674156 E 04

1.041759183084 E 04 1.057326037873 E 04 1.011790803983 E 04 1.049386278292 E 04 1.058378030620 E 04 1.054708986040 E 04 1.011790032703 E 04 1.108373250837 E 04 1.108480958884 E 04 1.022916214282 E 04 1.102338908514 E 04 1.108480958884 E 04 1.074600124223 E 04 3.794582738463 E 04 1.018629442582 E 04

1.256928753897 E 04 1.044372085311 E 04 1.264411600874 E 04 1.049701145029 E 04 1.044372085311 E 04 1.329069314639 E 04 1.002787557487 E 04 1.261988747883 E 04 1.261988747883 E 04 1.013595589892 E 04 1.114846685264 E 04 1.261988747883 E 04 1.265188041055 E 04 1.455053118626 E 04 1.480491820702 E 04

1.168400700815 E 04 1.070228566593 E 04 1.039802035425 E 04 1.115842157276 E 04 1.070749093697 E 04 1.825495131990 E 04 1.025352821628 E 04 1.660063460051 E 04 1.660063460051 E 04 1.441205094906 E 04 1.100527444996 E 04 1.660063460051 E 04 1.048114630651 E 04 2.703283204307 E 04 1.843542173805 E 04

1.119790612532 E 04 1.131713423785 E 04 1.109139622632 E 04 1.085673561603 E 04 1.005819331417 E 04 1.101988847359 E 04 1.006258723472 E 04 1.126177632509 E 04 1.208526345748 E 04 1.301656854635 E 04 1.044232500190 E 04 1.214066233230 E 04 1.223016607783 E 04 1.232044773762 E 04 1.220219100673 E 04

Probability values P_Break_EI

P_Break_EP

P_Break_WW

P_BreakWS

P_Trans_ESEP

P_Trans_EPES

1.065213990477 E 04 1.008138851590 E 04 1.007070772957 E 04 1.065389939474 E 04 1.008138851590 E 04 1.106047825272 E 04 1.010839020368 E 04 1.162275956024 E 04 1.162275956024 E 04 3.202382589825 E 04 1.351505215346 E 04 1.162275956024 E 04 1.163329266621 E 04 1.926030021959 E 04 2.320826650758 E 04

1.033092594998 E 04 1.009876865473 E 04 1.009459551782 E 04 1.268121367717 E 04 1.009878004556 E 04 1.193607308880 E 04 1.009878004556 E 04 1.774954402573 E 04 1.189006500254 E 04 1.815884941859 E 04 1.074043450829 E 04 1.774954402573 E 04 1.774954402573 E 04 1.274568498726 E 04 1.275174696018 E 04

1.201413262963 E 04 1.011184760200 E 04 1.092726240777 E 04 1.024667924657 E 04 1.745240138795 E 04 1.196300955805 E 04 1.021193775234 E 04 1.381201842186 E 04 1.539554325945 E 04 1.037816562659 E 04 3.480614106757 E 04 3.493840487408 E 04 8.346862553374 E 04 5.894278730730 E 04 8.310515465996 E 04

1.028419371369 E 04 1.006559131653 E 04 1.229055625028 E 04 1.015660754597 E 04 1.006436221293 E 04 1.025133185030 E 04 1.002854044511 E 04 1.004306774731 E 04 1.011136779011 E 04 1.449414921199 E 04 1.220757133175 E 04 1.011136779011 E 04 1.011136779011 E 04 1.133224981433 E 04 1.218297822105 E 04

1.322491332642E-03

1.003480002186 E 04 1.004083311881 E 04 1.003024455106 E 04 1.089252442911 E 04 1.006507382681 E 04 1.082292682391 E 04 1.105475934571 E 04 1.092666087111 E 04 1.099590308208 E 04 1.094192145349 E 04 1.096873241332 E 04 1.092666087111 E 04 1.107145242193 E 04 1.108153175440 E 04 1.075871909428 E 04

1.262574824539E-03 1.270238833297E-03 1.364062403227E-03 1.262574824539E-03 1.524677092450E-03 1.675600491021E-03 1.545863739501E-03 1.842659580112E-03 1.867651779550E-03 1.586592256990E-03 1.842086017785E-03 1.842086017785E-03 1.867417019172E-03 2.224934291931E-03

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

Table A6 Grid size: 310 310. Target % conversion of S to P: 23%. The mean is calculated over 100 runs for each of the parameter sets. Standard deviation is calculated over the mean. % S to P shows the percentage of mean S converted to P after completion of simulation. The probability values are shown correct up to 12 decimal places. This optimization run gave 10 optimized parameter sets. Mean [P]

Standard deviation

% S to P

779.92

83.61669158

20.289282

798.63

26.67221723

20.77601457

798.99

25.9893121

20.78537981

804.37

30.15001382

20.92533819

813.7

26.34982223

21.16805411

820.93

28.85045209

21.35613944

825.28

27.71539338

21.46930281

829.94

31.67091519

21.5905307

833.36

31.49395324

21.67950052

854.21

25.32627294

22.22190427

Mean [P] Standard deviation % S to P

Probability values P_Join_ES

P_Join_EI

P_Join_EP

P_Join_WW

P_Join_WS

P_Break_ES

1.025154414949 E 04 1.016879157672 E 04 1.024989553227 E 04 1.016726324813 E 04 1.025154414949 E 04 1.311799593681 E 04 1.025153478869 E 04 1.261966007452 E 04 1.011875941111 E 04 1.313799126141 E 04

1.047655465574 E 04 1.003258630305 E 04 1.004048784817 E 04 1.003052648823 E 04 1.008414405943 E 04 1.003359546734 E 04 1.045039800440 E 04 1.045843471348 E 04 1.003604519990 E 04 1.004034256149 E 04

1.007538198862 E 04 1.016042535954 E 04 1.013336970058 E 04 1.013048995549 E 04 1.021795768715 E 04 1.023591086511 E 04 1.021467327985 E 04 1.001532081263 E 04 1.054990020144 E 04 1.056375667270 E 04

1.029672900362 E 04 1.010467134400 E 04 1.008681627219 E 04 1.007854766346 E 04 1.034511575188 E 04 1.030190191181 E 04 1.001784019269 E 04 1.001784019269 E 04 1.007882984038 E 04 1.029635181947 E 04

1.784630521279 E 04 1.740404486406 E 04 1.323516566618 E 04 1.519333853128 E 04 1.633242494350 E 04 1.402102358539 E 04 1.772107291791 E 04 1.772107291791 E 04 1.191606185819 E 04 1.131475677046 E 04

1.009179961914 E 04 1.098010811110 E 04 1.009093222780 E 04 1.009093222780 E 04 1.009179961914 E 04 1.099625853093 E 04 1.009644900279 E 04 1.009644900279 E 04 1.014261064234 E 04 1.009179961914 E 04

Probability values P_Break_EI

779.92

83.61669158

798.63

26.67221723

798.99

25.9893121

804.37

30.15001382

813.7

26.34982223

820.93

28.85045209

825.28

27.71539338

829.94

31.67091519

833.36

31.49395324

854.21

25.32627294

20.289282

1.024219015342 E 04 20.77601457 1.013958579476 E 04 20.78537981 1.003832728477 E 04 20.92533819 1.003804383025 E 04 21.16805411 1.033325838548 E 04 21.35613944 1.030140194842 E 04 21.46930281 1.000474419996 E 04 21.5905307 1.001096451244 E 04 21.67950052 1.034713536692 E 04 22.22190427 1.001709166702 E 04

P_Break_EP

P_Break_WW

5.317831670002 E 04 5.298775349379 E 04 5.205791700196 E 04 5.318321571071 E 04 5.392583801367 E 04 7.473682592545 E 04 5.445237717103 E 04 5.634520153266 E 04 7.552987412185 E 04 7.560774706415 E 04

1.001330825296E-03

P_BreakWS

1.042089431762 E 04 9.960849085510 1.000395045312 E 04 E 04 7.577077816457 1.042038876732 E 04 E 04 9.126619490115 1.042042091788 E 04 E 04 9.971318077005 1.041810910203 E 04 E 04 5.256913834617 1.084946735508 E 04 E 04 1.003554793300E-03 1.005504308560 E 04 1.016067097805E-03 1.005504308560 E 04 5.418352991011 1.327612314742 E 04 E 04 7.595356954257 1.083656056692 E 04 E 04

P_Trans_ESEP

P_Trans_EPES

3.453746818991E-03

1.018048588433 E 04 3.520053841059E-03 1.015076505366 E 04 3.975525013760E-03 1.029002452220 E 04 3.975525013760E-03 1.221525902957 E 04 3.939718532383E-03 1.009253138523 E 04 4.941083299202E-03 1.034675162228 E 04 4.351382405809E-03 1.008709788831 E 04 4.351382405809E-03 1.016744172282 E 04 5.366328479781E-03 1.031047850422 E 04 5.403328504164E-03 1.018114162481 E 04

Table A7 Grid size: 110 110. Target % conversion of S to P: 90%. The mean is calculated over 100 runs for each of the parameter sets. Standard deviation is calculated over the mean. % S to P shows the percentage of mean S converted to P after completion of simulation. The probability values are shown correct up to 12 decimal places. This optimization run gave 5 optimized parameter sets. Mean [P]

408.8 412.31

Standard deviation

8.895578291 10.43275243

% S to P

84.46280992 85.18801653

422.22

8.721944923

87.23553719

423.8

9.2463479

87.56198347

427.44

9.491431543

88.31404959

Mean [P]

408.8 412.31

Standard deviation

8.895578291 10.43275243

% S to P

84.46280992 85.18801653

Probability values P_Join_ES

P_Join_EI

P_Join_EP

P_Join_WW

P_Join_WS

P_Break_ES

1.906366343960 E 04 3.793926680111 E 04 5.220558651861 E 04 5.220558651861 E 04 5.339205808110 E 04

2.168982915735 E 04 1.259738529049 E 04 1.737943900325 E 04 1.737943900325 E 04 1.737943900325 E 04

2.407758287537 E 04 2.407758287537 E 04 2.258342664676 E 04 2.258342664676 E 04 2.258032271810 E 04

7.890955295886 E 04 7.890955295886 E 04 7.171240388621 E 04 7.171240388621 E 04 7.044489960122 E 04

1.959335046228 E 04 2.443670516304 E 04 2.835382728411 E 04 2.835382728411 E 04 2.835382728411 E 04

4.574376782015 E 04 4.501279534938 E 04 1.495491678341 E 04 1.495491678341 E 04 1.498161283462 E 04

Probability values P_Break_EI

P_Break_EP

P_Break_WW

P_BreakWS

P_Trans_ESEP

P_Trans_EPES

1.346148906259 E 04 1.526464084332 E 04

2.809809920226 E 04 2.747373838246 E 04

1.109643117549 E 01 1.116030715742 E 01

3.930060568939 E 04 3.930060568939 E 04

5.748170814911 E 04 5.730942555388 E 04

5.962502879476 E 04 3.475146233411 E 04

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

13

Table A7 (continued ) Mean [P]

Standard deviation

% S to P

422.22

8.721944923

87.23553719

423.8

9.2463479

87.56198347

427.44

9.491431543

88.31404959

Probability values P_Break_EI

P_Break_EP

P_Break_WW

P_BreakWS

P_Trans_ESEP

P_Trans_EPES

3.278459628989 E 04 3.278459628989 E 04 3.278459628989 E 04

1.226821783219 E 04 1.226821783219 E 04 2.887858561739 E 04

1.685225256082 E 01 1.685225256082 E 01 1.733324553309 E 01

4.648766809812 E 04 4.648766809812 E 04 4.648766809812 E 04

3.820425536883 E 04 3.820425536883 E 04 3.820425536883 E 04

6.476953442893 E 04 6.476953442893 E 04 1.752056886860 E 04

Table A8 Grid size: 310 310. Target % conversion of S to P: 90%. The mean is calculated over 100 runs for each of the parameter sets. Standard deviation is calculated over the mean. % S to P shows the percentage of mean S converted to P after completion of simulation. The probability values are shown correct up to 12 decimal places. This optimization run gave 7 optimized parameter sets. Mean [P]

3312.49

Standard deviation

336.0249659

% S to P

86.17299688

3351.21

25.67642278

87.18028096

3372.23

26.97954893

87.72710718

3373.92

26.3736477

87.7710718

3390.75

26.62093617

88.20889698

3416.96

29.49607747

88.89073881

3419.86

27.34146959

88.96618106

Mean [P]

3312.49

Standard deviation

336.0249659

% S to P

86.17299688

3351.21

25.67642278

87.18028096

3372.23

26.97954893

87.72710718

3373.92

26.3736477

87.7710718

3390.75

26.62093617

88.20889698

3416.96

29.49607747

88.89073881

3419.86

27.34146959

88.96618106

Probability values P_Join_ES

P_Join_EI

P_Join_EP

P_Join_WW

P_Join_WS

P_Break_ES

2.782291840010 E 04 2.782291840010 E 04 4.180627499917 E 04 4.180627499917 E 04 4.681376467652 E 04 4.180627499917 E 04 2.587588943294 E 03

2.539182217565 E 04 2.539182217565 E 04 1.007544383397 E 04 1.835464909500 E 04 4.421274939147 E 04 1.007544383397 E 04 2.229012383606 E 04

2.668443872555 E 04 2.102698966003 E 03 2.278154607271 E 03 2.278154607271 E 03 2.155548389858 E 04 2.278154607271 E 03 2.243646577081 E 03

8.793389944909 E 02 8.793389944909 E 02 9.811537678993 E 02 9.811537678993 E 02 9.056598410979 E 02 9.811537678993 E 02 8.248789519626 E 02

1.458984853307 E 04 1.463185534894 E 04 1.363549873094 E 04 1.363549873094 E 04 1.014459279551 E 04 1.363549873094 E 04 1.284608885461 E 03

3.300332656750 E 04 3.326002456803 E 04 2.978323261390 E 04 3.309241297116 E 04 2.613542098276 E 04 2.978323261390 E 04 3.669342311070 E 04

P_Break_EI

P_Break_EP

P_Break_WW

P_BreakWS

P_Trans_ESEP

P_Trans_EPES

1.504071109364 E 04 1.504071109364 E 04 1.061100634339 E 03 2.190494023141 E 03 7.366455185358 E 04 1.061100634339 E 03 2.103390601504 E 03

1.361938669141 E 03 1.361938669141 E 03 1.536581784585 E 04 1.536581784585 E 04 5.995170265665 E 04 1.536581784585 E 04 1.338758909279 E 03

1.294014551757 E 01 1.294014551757 E 01 1.252822090203 E 01 1.252822090203 E 01 1.219578338835 E 01 1.252822090203 E 01 1.293671863198 E 01

4.699319427507 E 02 4.699319427507 E 02 4.759974962949 E 02 4.759974962949 E 02 7.387452007468 E 02 4.759974962949 E 02 4.733699909692 E 02

4.877093122919 E 02 4.898199401191 E 02 5.064988241132 E 02 5.001403395937 E 02 5.107045233533 E 02 5.064988241132 E 02 7.885583452665 E 02

1.085400903442 E 04 1.085400903442 E 04 1.058449881287 E 04 1.058449881287 E 04 1.191895756350 E 04 1.058449881287 E 04 1.049082369600 E 04

Probability values

the experimental work. For the pseudo-artiﬁcial reaction simulations, computational resources (Stevin Supercomputer Infrastructure) and services were used in this work and provided by Ghent University, the Hercules Foundation and the Flemish Government —Department EWI while for the model validation, computational resources of Amazon Web Service (AWS) were used to run parallel threaded simulations. The resources were available on credit basis with its elastic cloud infrastructure (EC2).

Appendix A See Tables A1–A8.

References Abkowitz, J., Catlin, S., Guttorp, P., 1996. Evidence that hematopoiesis may be a stochastic process in vivo. Nature Medicine 2, 190–197. Alemani, D., Pappalardo, F., Pennisi, A., Brusic, V., 2012. Combining cellular automata and lattice Boltzmann method to model multiscale avascular tumor growth coupled with nutrient diffusion and immune competition. Journal of Immunological Methods 376 (1–2), 55–68. Alvarez, E.F., 1961. The kinetic and mechanism of the hydrolysis of phosphoric acid esters by potato acid phosphatase. Biochimica et Biophysica Acta 59, 672–680. Apte, A., Cain, J.W., Bonchev, D.G., Fong, S.S., 2008. Cellular automata simulation of topological effects on the dynamics of feed-forward motifs. Journal of Biological Engineering 2 (2), 1–12. Blake, W.J., Kaern, M., Cantor, C.R., Collins, J.J., 2003. Noise in eukaryotic gene expression. Nature 422, 633–637. Bull, H., Murray, P.G., Thomas, D., Fraser, A.M., Nelson, P.N., 2002. Acid phosphatases. Molecular Pathology 55 (2), 65–72.

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

14

S. Kar et al. / Chemical Engineering Science ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Chopard, B., Masselot, A., 1999. Cellular automata and lattice Boltzmann methods: a new approach to computational ﬂuid dynamics and particle transport. Future Generation Computer Systems 16 (2–3), 249–257. Coello, C.A.C., 1999. A comprehensive survey of evolutionary-based multiobjective optimization techniques. Knowledge Information Systems 1 (3), 269–308. Coello, C.A.C., 2006. Evolutionary multi-objective optimization: a historical view of the ﬁeld. IEEE Computational Intelligence Magazine 1 (1), 28–36. Coello, C.A.C., Lamont, G.B., Veldhuizen, D.A.V., 2007. Evolutionary Algorithms for Solving Multi-Objective Problems. Genetic and Evolutionary Computation, second ed. Springer. Deb, K., Agrawal, R.B., 1995. Simulated binary crossover for continuous search space. Complex Systems 9 (2), 115–148. Deb, K., 2001. Multi Objective Optimization Using Evolutionary Algorithms. Wiley, New York. Deb, K., Pratap, A., Agarwal, A., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation 6, 182–197. Deb, K., Tiwari, S., 2006. Omni-optimizer: a generic evolutionary algorithm for single and multi-objective optimization. European Journal of Operational Research 185 (3), 1062–1087. Deutsch, A., Dormann, S., 2005. Cellular Automaton Modelling of Biological Pattern Formation. Birkhauser, Boston. Durillo, J.J., Nebro, A.J., 2011. jMetal: a Java framework for multi-objective optimization. Advances in Engineering Software 42 (10), 760–771. Dutta, A., Kar, S., Constales, D., Nopens, I., 2011. Modeling of ﬁrst-order enzyme kinetic reactions: a Cellular Automata approach. In: 7th International Workshop on Mathematics in Chemical Kinetics and Engineering (MaCKiE). Heidelberg, Germany. Dutta, A., Kar, S., Constales, D., Nopens, I. Modeling of ﬁrst-order enzyme kinetic reactions: a Cellular Automata approach, submitted for publication. Fagerlind, M.G., Webb, J.S., Barraud, N., McDougald, D., Jansson, A., Nilsson, P., Harlén, M., Kjelleberg, S., Rice, S.A., 2012. Dynamic modelling of cell death during bioﬁlm development. Journal of Theoretical Biology 295, 23–36. Goldberg, D.E., 1989. Genetic algorithms in search, Optimization and Machine Learning. Addison-Wesley Publishing Company, Reading, Massachusetts. Holland, J.H., 1975. Adaptation in Natural and Artiﬁcial Systems. MIT Press, Ann Arbor, MIT.

Kansal, A.R., Torquato, S., Harsh IV, G.R., Chiocco, E.A., Deisboeck, T.S., 2002. Simulated brain tumor growth dynamics using a three-dimensional cellular automaton. Journal of Theoretical Biology 203 (4), 367–382. Kar, S., Ganai, S., Dutta, A., Dutta, D., Chaudhuri, S., 2010. A sensitivity analysis study of enzyme inhibition kinetics through cellular automata. AIP conference Proceedings (ICMOC) 98 1, 301–306. Kier, L.B., Bonchev, D., Buck, G.A., 2005. Modeling biochemical networks: a cellularautomata approach. Chemistry and Biodiversity 2 (2), 233–243. Konak, A., Cott, D.W., Smith, A.E., 2006. Multi-objective Optimization Using Genetic Algorithms: A Tutorial, Reliability Engineering Systems Safety, vol. 91. Elsevier, pp. 992–1007. Lauf, S., Haase, S., Hostert, P., Lakes, T., Kleinschmit, B., 2012. Uncovering land-use dynamics driven by human decision-making—a combined model approach using cellular automata and system dynamics. Environmental Modeling and Software 27–28, 71–82. Meng, T.C., Somani, S., Dhar, P., 2004. Modeling and simulation of biological systems with stochasticity. In Silico Biology 4, 293–309. Nag, K., Pal, T., 2012. A new Archive-based Steady-state Genetic Algorithm. World Congress on Computational Intelligence, WCCI 2012 IEEE. Brisbane, Australia, pp. 853–859. Odds, F.C., Hierholzer, J.C., 1973. Puriﬁcation and properties of a glycoprotein acid phosphatase from Candida albicans. Journal of Bacteriology 114 (1), 257–266. Pizarro, G., Teixeira, J., Sepúlveda, M., Noguera, D., 2005. Bitwise implementation of a two-dimensional cellular automata bioﬁlm model. Journal of Computing in Civil Engineering 19 (3), 258–268. Somers, J.A., 1993. Direct simulation of ﬂuid ﬂow with cellular automata and the lattice-Boltzmann equation. Applied Scientiﬁc Research 51, 127–133. Soarse-Filho, B.S., Cerqueira, G.C., Pennachin, C.L., 2002. Dinamica—a stochastic cellular automata model designed to simulate the landscape dynamics in an Amazonian colonization frontier. Ecological Modeling 154 (3), 217–235. Wu, M., Higgs, P.G., 2012. The origin of life is a spatially localized stochastic transition. Biology Direct 24, 7–42. Zitzler, E., Laumanns, M., Thiele, L., 2001. SPEA2: improving the strength pareto evolutionary algorithm. Computational Engineering Network Laboratory (TIK). Swiss Federal Institute of Technology (ETH). Zurich, Switzerland. Technical Report 103.

Please cite this article as: Kar, S., et al., An improved cellular automata model of enzyme kinetics based on genetic algorithm. Chem. Eng. Sci. (2013), http://dx.doi.org/10.1016/j.ces.2013.08.013i

An improved cellular automata model of enzyme kinetics based on genetic algorithm

An improved cellular automata model of enzyme kinetics based on genetic algorithm

Recommend Documents