Team:Edinburgh/Cellulases (MATLAB model)

Cellulases (MATLAB model)

MATLAB is a general-purpose mathematical tool, developed by [http://www.mathworks.com Mathworks], and commonly used by engineers. Among other things, it can be used to solve ordinary differential equations (ODE's) by numerical integration.

An attempt was made to use MATLAB to model the degradation of cellulose into glucose in a biorefinery. But accurately predicting how much is converted in the synergistic system (where enzymes are kept close together) is difficult without literature to provide the ODE's and the kinetic parameters. Therefore this model only looks at the free floating (non-synergistic) enzyme approach. It is deterministic and is set by a series of initial conditions.

Assumptions

The mathematical model is based on the ODE's and kinetic parameters outlined in [http://onlinelibrary.wiley.com/doi/10.1021/bp034316x/full Kadam et al (2004)]. The following are its assumptions and basis:

• Underlying assumption - cellulose, cellobiose and glucose change continuously with time
• Rate equations assume enzyme adsorption follows the Langmuir isotherm model.
• Glucose and cellobiose, which are the products of cellulose hydrolysis, are assumed to "competitively inhibit enzyme hydrolysis".
• All reactions are assumed to follow the same temperature dependency Arrhenius relationship (shown below). However, in reality it should be different for every enzyme component, "because of their varying degrees of thermostability, with β-glucocidase being the most thermostable. Hence the assumption is a simplification of reality".
• Conversion of cellobiose to glucose follows the Michaelis-Menten enzyme kinetic model.

Equations

Rate Equations

Cellulose to cellobiose reaction with competitive glucose, cellobiose and xylose inhibition.

Cellulose to glucose reaction with competitive glucose, cellobiose and xylose inhibition.

Cellobiose to glucose reaction with competitive glucose, cellibiose and xylose inhibition.

Constants

k_nr — reaction rate constant for reaction n

E_nB — the bound concentration for exo and endo-β-1,4-glucanase for reaction n

R_s — substrate reactivity parameter

S — substrate reactivity at a given time (g/kg)  

G₂ — concentration of cellobiose

G — concentration of glucose  

X — xylose concentration  

K_nIG2 — inhibition constant for cellobiose at reaction n

K_nIG — inhibition constant for glucose at reaction n

K_nIX — xylose inhibition constant for reaction n

Note: For simplicity's sake we have assumed no xylose in the system, therefore X=0.

Langmuir Isotherm

The Langmuir Isotherm model mathematically describes enzyme adsorption onto solid cellulose substrates. Even though the Langmuir model is based on uniform binding sites and no interaction between the adsorbing molecules, it is not valid for cellulase adsorption onto cellulose. "Nevertheless the Langmuir formulation remains useful for mathematically describing the phenomenon of enzyme adsorption" ([http://onlinelibrary.wiley.com/doi/10.1021/bp034316x/full Kadam et al, 2004]).

Constants

exo and endo-β-1,4-glucanase, i=1  

β-glucosidase, i=2  

E_imax — Maximum mass of enzyme that can be absorbed onto a unit of mass substrate

K_iad — Dissociation constant for enzyme i

E_iF — Free enzyme concentration for enzyme i

S — Substrate reactivity at a given time (g/kg)

Mass Balances

Cellulose mass balance

Cellobiose mass balance

Glucose mass balance

Arrhenius Equation

Arrhenius equation is an empirical relationship which is used to model the temperature dependent reaction rate constant. Note: T1 is set at 45 °C

Constants

K_ir — Reaction rate constant of reaction i

E_ai — Activation energy of reaction i

R — Universal gas constant

Construction of Model

The model was constructed using the numerical programme MATLAB. A script file was generated which holds the variable dictionary, constants, temperature dependency equations, the "ODE45" differential equation solver, and the plot command.

A separate function file to the script is created as script files can only operate on the variables that are coded into their m-files. R_s (the substrate reactivity parameter) changes at every iteration because it is dependent on the S the substrate concentration at a given time, the substrate being cellulose. Therefore S at the first iteration is S₀, the initial substrate concentration. At the second iteration S is the previous value calculated by the ODE, and so on. After each step the new value of S is fed into the function file and is used for calculating the reaction rate constant for cellobiose and glucose, etc.

R_s — substrate reactivity equation

ODE45 calls on the function file to calculate the rate equations and then substitute them into the respective mass balance. A numerical integration is performed and the results can be seen below. ODE45 is used as it is more accurate than other solvers [http://www.mathworks.co.uk/help/techdoc/ref/ode23.html (according to Mathworks)]. It is based on an explicit Runge-Kutta formula and is a one-step solver.

Results

Figure 1: Graph of cellulose degradation over time with maximum β-glucosidase.
Initial conditions: Cellulose - 100 g/kg, Glucose - 0.01 g/kg , Cellobiose - 0.01 g/kg.
Enzymes: Exo/endo-glucanase - 0.01 g/kg, β-glucosidase - 1 g/kg
Temperature 35° C

Figure 2: Graph of cellulose degradation over time with maximum Exo/endo-glucanase.
Initial conditions: Cellulose - 100 g/kg, Glucose - 0.01 g/kg , Cellobiose - 0.01 g/kg.
Enzymes: Exo/endo-glucanase - 1 g/kg, β-glucosidase - 0.01 g/kg
Temperature 35° C

Figure 1 is set with β-glucosidase at its maximum concentration and Figure 2 with Exo/endo-glucanase at maximum. This is to compare the effect of certain enzymes on cellulose degradation and glucose production. The result is consistent with what is expected. Exoglucanase chews away at the end of a cellulose chain, producing cellobiose sugars and endoglucanase cuts cellulose chanins in the centre, turning one chain into two. This concurs with the results, with Figure 2 modelling markedly higher amount of cellulose at the maximum amount of exo/endo-glucanase than that of Figure 1. Whereas β-glucosidase cuts cellobiose in half, producing two glucose molecules. Figure 1 which β-glucosidase is at its maximum produces 167% more glucose than in Figure 2 with β-glucosidase at its minimum at steady state.

Cellobiose in Figure 1 increases to 10 g/kg in 20 hours but after 50 hours begins to decrease. This is because β-glucosidase starts to convert it to glucose. This result cannot be seen in Figure 2 with little β-glucosidase present in the system as a result cellobiose increases to 60 g/kg.

Figure 3: Graph of cellulose degradation over time with maximum β-glucosidase.
Initial conditions: Cellulose - 100 g/kg, Glucose - 0.01 g/kg , Cellobiose - 0.01 g/kg.
Enzymes: Exo/endo-glucanase - 0.01 g/kg, β-glucosidase - 1 g/kg
Temperature 35° C. The x-axis is set to semilog.

Figure 3 is under the same conditions as Figure 1, but its x-axis is set to a semilog scale to illustrate what happens over a longer period of time. The effect of high β-glucosidase can better be seen with the complete conversion of cellobiose to glucose.

Interestingly, the amount of glucose over time highlights one of the limitations of an ODE based model. The amount of glucose increases beyond the amount of cellulose present initially, which breaks the law of conservation of mass. This is caused by continuous iterations by MATLAB which don't stop after a certain point in time. The ODE therefore has limits within itself where it works best. This is one reason why Edinburgh decided to construct an agent-based model, using Kappa. That model explicitly tracks individual sugar molecules and so avoids this problem.

Download MATLAB file

You can download the ready model as a .zip file:

References

• Kadam KL, Rydholm EC, McMillan JD (2004) [http://onlinelibrary.wiley.com/doi/10.1021/bp034316x/full Development and Validation of a Kinetic Model for Enzymatic Saccharification of Lignocellulosic Biomass]. Biotechnology Progress 20(3): 698–705 (doi: 10.1021/bp034316x).