ORIGINAL PAPER

Occurrence of dead core in catalytic particles containing
immobilized enzymes: analysis for the Michaelis–Menten kinetics
and assessment of numerical methods

Félix Monteiro Pereira1 • Samuel Conceição Oliveira2

Received: 26 March 2016 / Accepted: 23 June 2016 / Published online: 30 June 2016

� Springer-Verlag Berlin Heidelberg 2016

Abstract In this article, the occurrence of dead core in

catalytic particles containing immobilized enzymes is

analyzed for the Michaelis–Menten kinetics. An assess-

ment of numerical methods is performed to solve the

boundary value problem generated by the mathematical

modeling of diffusion and reaction processes under steady

state and isothermal conditions. Two classes of numerical

methods were employed: shooting and collocation. The

shooting method used the ode function from Scilab soft-

ware. The collocation methods included: that implemented

by the bvode function of Scilab, the orthogonal collocation,

and the orthogonal collocation on finite elements. The

methods were validated for simplified forms of the

Michaelis–Menten equation (zero-order and first-order

kinetics), for which analytical solutions are available.

Among the methods covered in this article, the orthogonal

collocation on finite elements proved to be the most robust

and efficient method to solve the boundary value problem

concerning Michaelis–Menten kinetics. For this enzyme

kinetics, it was found that the dead core can occur when

verified certain conditions of diffusion–reaction within the

catalytic particle. The application of the concepts and

methods presented in this study will allow for a more

generalized analysis and more accurate designs of hetero-

geneous enzymatic reactors.

Keywords Dead core � Michaelis–Menten kinetics �
Diffusion and reaction � Immobilized enzymes �
Effectiveness factor

Introduction

The main objectives of the development of new processes

employing biological agents, such as cells and enzymes,

include the manufacturing of new products, improvements

in product quality, the minimization of environmental

impacts, and cost reductions [1, 2].

The bioprocesses include the transformation of raw

materials (substrates) in products into bioreactors using

cells or enzymes. The selection of the transformation agent

and of the bioreactor operation mode depends on the

assessment of the advantages and disadvantages presented

by all possible configurations [1, 2].

The bioreactors containing immobilized enzymes are

employed in several industrial processes, such as: the

treatment of wastewater and the production of pharma-

ceuticals, chemicals, foods, beverages, biofuels, enzymes,

among other bioproducts [1, 2].

The application of free enzymes within a bioreactor can

cause a biocatalyst loss, which exits the reactor through the

outlet flow. One of the limiting factors when using

enzymes at an industrial scale is the cost of these enzymes,

which can be reduced if the enzymes can be reused. The

most feasible possibility of reuse is found in the immobi-

lization of the enzymes in matrices [3].

The advantages resulting from enzyme immobilization

have led to the development of new industrial processes

& Samuel Conceição Oliveira

samueloliveira@fcfar.unesp.br

Félix Monteiro Pereira

felix@dequi.eel.usp.br

1 Departamento de Engenharia Quı́mica, Escola de Engenharia

de Lorena, USP–Universidade de São Paulo, Estrada

Municipal do Campinho, 12602-810 Lorena, SP, Brazil

2 Departamento de Bioprocessos e Biotecnologia, Faculdade

de Ciências Farmacêuticas, UNESP–Univ Estadual Paulista,

Rodovia Araraquara-Jaú Km 1, 14800-903 Araraquara, SP,

Brazil

123

Bioprocess Biosyst Eng (2016) 39:1717–1727

DOI 10.1007/s00449-016-1647-0

http://crossmark.crossref.org/dialog/?doi=10.1007/s00449-016-1647-0&domain=pdf
http://crossmark.crossref.org/dialog/?doi=10.1007/s00449-016-1647-0&domain=pdf


employing heterogeneous enzymatic reactors. The model-

ing and simulation of these bioreactors is one of main steps

for the design, optimization, and control of the entire

process of conversion. In most cases, when modeling and

simulating enzymatic reactors, it is assumed that the

operating conditions are ideal, that is: plug or mixed flow

pattern, constant temperature, and constant pH. Further-

more, partitioning effects, diffusion resistance and thermal

inactivation of the biocatalyst are commonly neglected in

conventional modeling. However, in some cases, immo-

bilization can cause a high diffusion resistance, requiring a

detailed description of the diffusion–reaction processes

within the biocatalyst.

When the diffusion resistance is quite high, the substrate

concentration can be equal to zero in a given biocatalyst

region, as shown in Fig. 1. In this region, known as the

dead core, the reaction does not occur. Thus, the designed

reactors that assume that all catalyst volume is active will

present a lower performance than what is expected, justi-

fying, in such cases, the use of an appropriate mathematical

model for an accurate design of these reactors, which are

subjected to mass transfer limitations.

Analysis of the dead-core occurrence for zero- and first-

order kinetics can be carried out analytically. For the

Michaelis–Menten kinetics, due to its non linearity, the

analysis can only be conducted by applying appropriate

numerical methods. The Michaelis–Menten equation suc-

cessfully describes the kinetic behavior of a large number

of reactions catalyzed by enzymes, thus justifying its study.

Aiming to contribute through the design and operation

of bioreactors with immobilized enzymes, this study pre-

sents a methodology to analyze the dead-core occurrence

for an enzymatic reaction whose kinetics is described by

the Michaelis–Menten equation. For this aim, classical

numerical methods were tested for the solution of boundary

value problem generated by the theoretical analysis of

diffusion and reaction processes.

To test the proposed numerical methodology, simplified

forms of Michaelis–Menten kinetics, corresponding to

zero- and first-order kinetics, were analyzed. For these

cases, the analytical solutions are known for classical

geometries of biocatalysts (slab, cylinder, and sphere), thus

allowing one to draw comparisons with the numerical

solutions. After having validated the numerical methods,

the reaction–diffusion problem associated with Michaelis–

Menten kinetics was solved.

Mathematical modeling

The mass balance of substrate inside the particle, for steady

state and isothermal conditions, is given by [4–6]:

Def

d2S

dX2
þ a� 1ð Þ

X

dS

dX

� �
¼ v ð1Þ

In Eq. 1, S is the substrate concentration, X is the

spatial coordinate, Def is the effective-diffusion coeffi-

cient of substrate inside the particle, v is the rate of

substrate consumption, and a is the geometric shape

factor of the particle (a = 1 for slab, a = 2 for cylinder

and a = 3 for sphere).

Equation 1 can be given in terms of dimensionless

parameters and variables, as follows:

d2s

dx2
þ a� 1ð Þ

x

ds

dx
¼ a2/2 vðSÞ

vðS0Þ
ð2Þ

where

x ¼ X

XL

ð3Þ

s ¼ S

S0

ð4Þ

/2 ¼ X2
L

a2

v S0ð Þ
DefS0

ð5Þ

In Eqs. 2, 3, 4 and 5, x is the dimensionless spatial

coordinate, s is the dimensionless substrate concentration,

XL is the half value of thickness for slab and the radius for
Fig. 1 Representation of dead core in a spherical particle of a

biocatalyst

1718 Bioprocess Biosyst Eng (2016) 39:1717–1727

123



cylinder and sphere, S0 is the substrate concentration on the

particle surface, and / is the Thiele modulus.

Equation (2) is subjected to the following boundary

conditions:

ds

dx
¼ 0 at x

¼ a; for /[/crit problem with dead-core occurrenceð Þ
ð6Þ

ds

dx
¼ 0 at x

¼ 0; for /�/crit problem without dead-core occurrenceð Þ
ð7Þ

s ¼ 1 at x ¼ 1 ð8Þ

In Eq. 6, a is the dead-core position (dimensionless),

and /crit is the critical value of the Thiele modulus above

which the dead core occurs.

Equation 9 can be used to calculate the internal effec-

tiveness factor (g) [6, 7], a parameter ranging from 0 to 1

that indicates how much diffusion resistance controls the

rate of reaction.

g ¼ 1

a/2

ds

dx

����
x¼1

ð9Þ

Equation (10) shows the substrate consumption rate

given by Michaelis–Menten kinetics [4, 5]:

v Sð Þ ¼ vmaxS

K þ S
ð10Þ

where, K is the Michaellis–Menten constant, and vmax is the

maximum rate of substrate consumption.

The Michaelis–Menten equation has the following

asymptotic properties [4, 5]:

(a) For low substrate concentrations (s � K), it can be

approximated to first-order kinetics:

v Sð Þ ¼ k1S ; where k1 ¼ vmax=K ð11Þ

(b) For high substrate concentrations (s � K), it can be

approximated to zero-order kinetics:

v Sð Þ ¼ k0 ; where k0 ¼ vmax ð12Þ

Thus, introducing Eqs. 10, 11 and 12 into Eq. (2):

Fig. 2 Flowchart to calculate the concentration profiles, effectiveness

factor, and dead-core position

Fig. 3 Flowchart for the numerical determination of the critical

Thiele modulus (/crit)

Bioprocess Biosyst Eng (2016) 39:1717–1727 1719

123



d2s

dx2
þ a� 1ð Þ

x

ds

dx
¼ a2/2 1 þ bð Þs

bþ s
ð13Þ

d2s

dx2
þ a� 1ð Þ

x

ds

dx
¼ a2/2s ð14Þ

d2s

dx2
þ a� 1ð Þ

x

ds

dx
¼ a2/2 ð15Þ

In Eq. 13, b is the dimensionless Michaelis–Menten

constant (b = K/S0).

The exact analytical solutions of Eqs. 14 and 15 can be

found using the software Mathematica [8]. These solutions

were used to assess the accuracy of the numerical methods

applied in this study to solve Eq. 13.

The analytical solutions for the zero-order kinetics

(Eq. 15), subjected to the boundary conditions given by

Eqs. 6, 7 and 8, for slab (Eqs. 16 and 17), cylinder (Eqs. 18

and 19), and sphere (Eqs. 20 and 21) particles, are pre-

sented as follows.

s ¼ 1

2
2 þ x2 � 1 � 2a x � 1ð Þ

� �
/2

� �
; for /[/crit

ð16Þ

s ¼ 1

2
2 þ x2 � 1

	 

/2

� �
; for /�/crit ð17Þ

s ¼ 1 þ x2 � 1
	 


/2 � 2a2/2 ln xð Þ; for /[/crit ð18Þ

s ¼ 1 þ x2 � 1
	 


/2; for /�/crit ð19Þ

Table 1 Results from the

shooting method for the BVP

solution for zero- and first-order

kinetics

Order a / /a /n error/ a ga gn errorg errors

0 1 0.5
ffiffiffi
2

p
1.4142 1.4 9 10-13 0 1 1 0 2.5 9 10-14

0 1 1
ffiffiffi
2

p
1.4142 1.4 9 10-13 0 1 1 0 9.9 9 10-14

0 1 2
ffiffiffi
2

p
1.4142 1.4 9 10-13 0.2929 0.7071 0.7071 8.8 9 10-13 2.5 9 10-12

0 1 4
ffiffiffi
2

p
1.4142 1.4 9 10-13 0.6464 0.3536 0.3536 1.1 9 10-13 6.2 9 10-13

0 1 8
ffiffiffi
2

p
1.4142 1.4 9 10-13 0.8232 0.1768 0.1768 1.4 9 10-14 1.6 9 10-13

0 2 0,5 1 1.0000 0 0 1 1 0 5.0 9 10-14

0 2 1 1 1.0000 0 0 1 1 1.1 9 10-10 9.0 9 10-11

0 2 2 1 1.0000 0 0.6184 0.6176 0.6176 4.8 9 10-9 4.3 9 10-9

0 2 4 1 1.0000 0 0.8173 0.3320 0.3320 6.3 9 10-10 1.0 9 10-8

0 2 8 1 1.0000 0 0.9102 0.1715 0.1715 6.4 9 10-10 8.3 9 10-9

0 3 0.5
ffiffiffi
6

p �
3 0.8165 8.1 9 10-14 0 1 1 3.7 9 10-14 7.7 9 10-14

0 3 1
ffiffiffi
6

p �
3 0.8165 8.1 9 10-14 0.3870 0.9421 0.9421 5.0 9 10-8 1.1 9 10-8

0 3 2
ffiffiffi
6

p �
3 0.8165 8.1 9 10-14 0.7409 0.5934 0.5934 1.4 9 10-8 8.3 9 10-9

0 3 4
ffiffiffi
6

p �
3 0.8165 8.1 9 10-14 0.8770 0.3255 0.3255 4.8 9 10-10 8.9 9 10-9

0 3 8
ffiffiffi
6

p �
3 0.8165 8.1 9 10-14 0.9399 0.1698 0.1698 2.7 9 10-10 5.4 9 10-9

1 1 0.5 a b a 0 0.9242 0.9242 8.4 9 10-8 9.5 9 10-8

1 1 1 a b a 0 0.7616 0.7616 5.0 9 10-9 7.1 9 10-8

1 1 2 a b a 0 0.4820 0.4820 5.9 9 10-9 5.3 9 10-8

1 1 4 a b a 0 0.2498 0.2498 4.9 9 10-11 1.5 9 10-8

1 1 8 a b a 0 0.1250 0.1250 2.8 9 10-14 2.7 9 10-8

1 2 0.5 a b a 0 0.8928 0.8928 5.9 9 10-9 3.6 9 10-8

1 2 1 a b a 0 0.6978 0.6978 1.3 9 10-8 1.1 9 10-7

1 2 2 a b a 0 0.4318 0.4318 9.0 9 10-10 2.2 9 10-8

1 2 4 a b a 0 0.2338 0.2338 2.8 9 10-10 3.1 9 10-8

1 2 8 a b a 0 0.1210 0.1210 3.6 9 10-11 1.8 9 10-8

1 3 0.5 a b a 0 0.8762 0.8762 1.4 9 10-8 4.1 9 10-8

1 3 1 a b a 0 0.6716 0.6716 2.6 9 10-8 1.8 9 10-7

1 3 2 a b a 0 0.4167 0.4167 9.8 9 10-9 2.1 9 10-7

1 3 4 a b a 0 0.2292 0.2292 1.4 9 10-9 1.6 9 10-7

1 3 8 a b a 0 0.1198 0.1198 4.3 9 10-11 3.0 9 10-8

a Dead-core does not occur for the analytical solution of first-order kinetics [7]
b Method is divergent for high values of Thiele’s modulus (determination of /n is not possible)

1720 Bioprocess Biosyst Eng (2016) 39:1717–1727

123



s ¼ 1 �
3 2a3 x� 1ð Þ þ x� x3
� �

/2

2x
; for /[/crit ð20Þ

s ¼ 1 � 3

2
1 � x2
	 


/2; for /�/crit ð21Þ

Since the dead core occurs for /[/crit, the values of

/crit for the zero-order kinetics determined analytically are

as follows: /crit ¼
ffiffiffi
2

p
(slab), /crit = 1 (cylinder), and

/crit ¼
ffiffiffi
6

p �
3 (sphere). The value of a is determined by

taking s = 0 and x = a in Eqs. 16, 18 and 20.

For zero-order kinetics, the effectiveness factor can be

calculated analytically by Eq. 22 for the three classical

geometries.

g ¼ 1 � aa ð22Þ

The analytical solutions for first-order kinetics (Eq. 14),

subjected to the boundary conditions given by Eqs. 7 and

8, are presented as follows for slab (Eq. 23), cylinder

(Eq. 24), and sphere (Eq. 25) particles [8–10].

s ¼ cosh x/ð Þ
coshð/Þ ð23Þ

s ¼ I0 2x/ð Þ
I0 2/ð Þ ð24Þ

s ¼ sinhð3x/Þ
x sinh 3/ð Þ ð25Þ

For first-order kinetics, the effectiveness factor can be

calculated analytically as follows for slab (Eq. 26), cylin-

der (Eq. 27), and sphere (Eq. 28) [8–10].

Table 2 Results from the

orthogonal collocation method

for the BVP solution for zero-

order and first-order kinetics

Order a / /a /n error/ a ga gn errorg errors

0 1 0.5
ffiffiffi
2

p
1.4142 2.9 9 10-8 0 1 1 2.5 9 10-6 6.0 9 10-7

0 1 1
ffiffiffi
2

p
1.4142 2.9 9 10-8 0 1 1 5.0 9 10-7 4.0 9 10-7

0 1 2
ffiffiffi
2

p
1.4142 2.9 9 10-8 0.2929 0.7071 0.7071 2.0 9 10-7 3.5 9 10-8

0 1 4
ffiffiffi
2

p
1.4142 2.9 9 10-8 0.6464 0.3536 0.3536 8.4 9 10-8 3.5 9 10-8

0 1 8
ffiffiffi
2

p
1.4142 2.9 9 10-8 0.8232 0.1768 0.1768 4.2 9 10-8 3.5 9 10-8

0 2 0.5 1 1.0000 1.0 9 10-7 0 1 1 8.4 9 10-6 3.0 9 10-7

0 2 1 1 1.0000 1.0 9 10-7 0 1 1 2.0 9 10-6 2.0 9 10-7

0 2 2 1 1.0000 1.0 9 10-7 0.6184 0.6176 0.6176 2.2 9 10-5 7.5 9 10-5

0 2 4 1 1.0000 1.0 9 10-7 0.8173 0.3320 0.3320 3.9 9 10-6 2.6 9 10-5

0 2 8 1 1.0000 1.0 9 10-7 0.9102 0.1715 0.1715 7.0 9 10-7 1.1 9 10-5

0 3 0.5
ffiffiffi
6

p �
3 0.8165 4.0 9 10-8 0 1 1 4.4 9 10-6 2.0 9 10-7

0 3 1
ffiffiffi
6

p �
3 0.8165 4.0 9 10-8 0.3814 0.9421 0.9418 2.2 9 10-4 8.1 9 10-4

0 3 2
ffiffiffi
6

p �
3 0.8165 4.0 9 10-8 0.7402 0.5934 0.5933 4.5 9 10-5 1.7 9 10-4

0 3 4
ffiffiffi
6

p �
3 0.8165 4.0 9 10-8 0.8769 0.3255 0.3255 1.0 9 10-5 6.4 9 10-5

0 3 8
ffiffiffi
6

p �
3 0.8165 4.0 9 10-8 0.9398 0.1698 0.1698 2.7 9 10-6 2.8 9 10-5

1 1 0.5 a b a 0 0.9242 0.9242 2.8 9 10-6 5.0 9 10-7

1 1 1 a b a 0 0.7616 0.7616 7.0 9 10-7 4.0 9 10-7

1 1 2 a b a 0 0.4820 0.4820 1.0 9 10-7 5.0 9 10-6

1 1 4 a b a 0 0.2498 0.2498 5.4 9 10-9 2.5 9 10-5

1 1 8 a b a 0 0.1250 0.1250 1.5 9 10-8 1.3 9 10-4

1 2 0.5 a b a 0 0.8928 0.8928 8.2 9 10-6 5.0 9 10-7

1 2 1 a b a 0 0.6978 0.6978 2.0 9 10-6 2.8 9 10-6

1 2 2 a b a 0 0.4318 0.4318 5.0 9 10-7 1.7 9 10-5

1 2 4 a b a 0 0.2338 0.2338 1.0 9 10-7 1.0 9 10-4

1 2 8 a b a 0 0.1210 0.1210 4.0 9 10-7 4.8 9 10-4

1 3 0.5 a b a 0 0.8762 0.8762 4.2 9 10-6 6.0 9 10-7

1 3 1 a b a 0 0.6716 0.6716 1.2 9 10-6 5.2 9 10-6

1 3 2 a b a 0 0.4167 0.4167 4.0 9 10-7 3.8 9 10-5

1 3 4 a b a 0 0.2292 0.2292 2.0 9 10-7 2.1 9 10-4

1 3 8 a b a 0 0.1198 0.1198 1.4 9 10-5 9.9 9 10-4

a Dead-core does not occur for the analytical solution of first-order kinetics [7]
b Converges to the maximum value of /n (/n = 100) for the bisection method

Bioprocess Biosyst Eng (2016) 39:1717–1727 1721

123



g ¼ tanh /ð Þ
/

ð26Þ

g ¼ I1ð2/Þ
/ I0ð2/Þ

ð27Þ

g ¼ 1

/ tanh 3/ð Þ �
1

3/2
ð28Þ

In Eqs. 24 and 27, I0 and I1 are the modified Bessel

functions of first kind and zero (I0) and one (I1) orders.

For the dead-core occurrence problem, incorporating the

Michaelis–Menten kinetics, numerical and exact analytical

solutions are not reported in the literature. For problems

without dead-core occurrence, some approximated analyt-

ical solutions are proposed [11–14]. As the dead core can

occur for high values of the Thiele modulus, a reliable

numerical method must be used for both the calculus of

internal effectiveness and the search/validation of approx-

imated analytical solutions for the problems including

dead-core occurrence.

In this article, different numerical methods were applied

to calculate the substrate concentration profile and internal

effectiveness factor for Michaelis–Menten kinetics, taking

the dead-core occurrence into account.

Numerical methodology

A computer program was developed in the Scilab software

to solve the boundary value problem (BVP) given by

Eqs. 2, 6, 7 and 8, and several numerical methods,

including shooting and collocation methods, were tested.

The shooting methods used the fsolve function of Scilab,

which is based on Newton’s method, coupled with the ode

function, which is the standard function to integrate

explicit ordinary differential equation systems. The ode

function is an interface to several solvers. Hence, the pre-

sent study tested the following methods:

• Nonstiff predictor–corrector Adams [15–17];

• Stiff backward differentiation formula [15–17];

• Adaptive Runge–Kutta of 4th order [18–20];

• Shampine and Watts program based on Fehlberg’s

Runge–Kutta pair of orders 4 and 5 [18–20].

Among the collocation methods used in this study, that

implemented by the bvode function of Scilab solves a

multi-point boundary value problem for a mixed order

system of ordinary differential equations, using algorithms

from the literature [21–24]. Other collocation methods used

in this study included the methods of orthogonal colloca-

tion and orthogonal collocation on finite elements [25, 26],

which were implemented in Scilab using the following

functions:

• legendre: to obtain the Legendre polynomials;

• lsqrsolve, poly, and roots: to calculate the collocation

points;

• fsolve (Newton’s method): to calculate the substrate

concentration on collocation points.

Figure 2 shows the flowchart to calculate the concen-

tration profile, effectiveness factor, and dead-core position.

Figure 3 presents the flowchart for the numerical

determination of the critical Thiele modulus.

The BVP solution method coupled with the bisection

method was applied to determine the dead-core position

(a) and the critical Thiele modulus (/crit).

In Figs. 2 and 3, aL, aR, /L, and /R are the minimum or

left (L-subscript) and maximum or right (R-subscript)

boundary values used in the bisection method to determine

a and /crit parameters.

Based on the flowcharts presented in Figs. 2 and 3, three

numerical methods were applied to the solution of the

BVP: shooting, orthogonal collocation, and orthogonal

collocation on finite elements.

Results and discussion

From the results obtained using the developed computer

programs, the numerical and analytical solutions for the

simplified forms of Michaelis–Menten equation were

compared to assess the validity and accuracy of the

Fig. 4 Profile of substrate concentration (/ = 50, a = 3, first-order

kinetics) calculated by orthogonal collocation method and orthogonal

collocation method on finite elements

1722 Bioprocess Biosyst Eng (2016) 39:1717–1727

123



proposed numerical methodologies. For this comparison,

the following parameters were analyzed:

• Internal effectiveness factor calculated by numerical

and analytical solutions and the absolute difference

between them: gn, ga, errorg = |ga - gn|;
• Critical Thiele modulus calculated by numerical and

analytical solutions and the absolute difference between

them: /n, /a, error/ = |/a - /n|;

• Maximum absolute difference between the values of

dimensionless substrate concentration calculated by

numerical and analytical solutions: errors

The results for each applied method are presented sep-

arately as follows.

Shooting methods

Table 1 shows the results obtained from the shooting

methods. The results were the same for all available

methods in the ode function of Scilab.

Results on Table 1 indicate that the shooting meth-

ods coupled with the ode function of Scilab were able

to solve the BVP for a wide range of Thiele modulus.

However, for high values of Thiele modulus, all

methods available in the ode function proved to be

divergent when the first-order kinetics was analyzed.

This result suggests that another methodology must be

employed to solve the nonlinear boundary value

problems.

Table 3 Results from

orthogonal collocation method

on finite elements for the BVP

solution for zero- and first-order

kinetics

Order a / /a /n error/ a ga gn errorg errors

0 1 0.5
ffiffiffi
2

p
1.4142 2.2 9 10-11 0 1 1 5.7 9 10-11 9.9 9 10-12

0 1 1
ffiffiffi
2

p
1.4142 2.2 9 10-11 0 1 1 1.1 9 10-11 6.9 9 10-12

0 1 2
ffiffiffi
2

p
1.4142 2.2 9 10-11 0.2929 0.7071 0.7071 2.6 9 10-11 2.7 9 10-12

0 1 4
ffiffiffi
2

p
1.4142 2.2 9 10-11 0.6464 0.3536 0.3536 2.7 9 10-11 2.6 9 10-12

0 1 8
ffiffiffi
2

p
1.4142 2.2 9 10-11 0.8232 0.1768 0.1768 2.8 9 10-11 2.8 9 10-12

0 2 0.5 1 1.0000 1.3 9 10-11 0 1 1 9.0 9 10-12 4.0 9 10-12

0 2 1 1 1.0000 1.3 9 10-11 0 1 1 2.0 9 10-12 2.0 9 10-12

0 2 2 1 1.0000 1.3 9 10-11 0.6184 0.6176 0.6176 3.6 9 10-11 2.8 9 10-12

0 2 4 1 1.0000 1.3 9 10-11 0.8173 0.3320 0.3320 3.3 9 10-11 2.6 9 10-12

0 2 8 1 1.0000 1.3 9 10-11 0.9102 0.1715 0.1715 3.8 9 10-11 2.7 9 10-12

0 3 0.5
ffiffiffi
6

p �
3 0.8165 2.7 9 10-12 0 1 1 2.9 9 10-12 2.6 9 10-12

0 3 1
ffiffiffi
6

p �
3 0.8165 2.7 9 10-12 0.3869 0.9421 0.9421 8.8 9 10-14 7.1 9 10-11

0 3 2
ffiffiffi
6

p �
3 0.8165 2.7 9 10-12 0.7402 0.5934 0.5934 1.1 9 10-12 2.6 9 10-12

0 3 4
ffiffiffi
6

p �
3 0.8165 2.7 9 10-12 0.8769 0.3255 0.3255 7.7 9 10-12 2.6 9 10-12

0 3 8
ffiffiffi
6

p �
3 0.8165 2.7 9 10-12 0.9398 0.1698 0.1698 3.0 9 10-11 2.6 9 10-12

1 1 0.5 a 86.9590 a 0 0.9242 0.9242 5.1 9 10-11 9.1 9 10-12

1 1 1 a 86.9590 a 0 0.7616 0.7616 8.4 9 10-12 5.7 9 10-12

1 1 2 a 86.9590 a 0 0.4820 0.4820 8.5 9 10-13 1.7 9 10-12

1 1 4 a 86.9590 a 0 0.2498 0.2498 7.6 9 10-14 3.0 9 10-11

1 1 8 a 86.9590 a 0 0.1250 0.1250 1.1 9 10-12 3.3 9 10-9

1 2 0.5 a 43.9220 a 0 0.8928 0.8928 7.0 9 10-12 3.4 9 10-12

1 2 1 a 43.9220 a 0 0.6978 0.6978 1.1 9 10-12 1.7 9 10-12

1 2 2 a 43.9220 a 0 0.4318 0.4318 1.2 9 10-13 1.8 9 10-11

1 2 4 a 43.9220 a 0 0.2338 0.2338 2.1 9 10-12 2.4 9 10-9

1 2 8 a 43.9220 a 0 0.1210 0.1210 8.5 9 10-10 2.6 9 10-7

1 3 0.5 a 29.5413 a 0 0.8762 0.8762 2.0 9 10-12 2.2 9 10-12

1 3 1 a 29.5413 a 0 0.6716 0.6716 2.3 9 10-13 1.4 9 10-12

1 3 2 a 29.5413 a 0 0.4167 0.4167 2.4 9 10-13 2.2 9 10-10

1 3 4 a 29.5413 a 0 0.2292 0.2292 1.0 9 10-10 3.1 9 10-8

1 3 8 a 29.5413 a 0 0.1198 0.1198 3.3 9 10-8 3.1 9 10-6

a Dead-core does not occur for the analytical solution of the first-order kinetics [7]

Bioprocess Biosyst Eng (2016) 39:1717–1727 1723

123



Collocation methods

Similar results (Table 2) were obtained using two collo-

cation methods: the orthogonal collocation and that

implemented by the bvode function of Scilab.

Table 2 shows a low accuracy to define the internal

effectiveness factor at high values of the Thiele modulus.

For high values from the Thiele modulus, the calculated

values of effectiveness factor do not agree with the ana-

lytical values. This fact was associated with the inconsis-

tent concentration profiles obtained for these cases, as

clearly shown in Fig. 4 regarding, orthogonal collocation.

Similar profiles were obtained for both the orthogonal

collocation method and that implemented by the bvode

function of Scilab.

Results presented in Table 2 and Fig. 4 regarding

orthogonal collocation suggest that another methodology

must be employed to solve nonlinear BVP.

The results for orthogonal collocation methods on finite

elements are presented in Table 3.

Table 3 shows an excellent accuracy for the calculated

values of effectiveness and substrate concentration. The

numerical values for the critical Thiele modulus also pre-

sent good accuracy when analyzing zero-order kinetics.

These calculations used 10 internal collocation points

(orthogonal collocation method) within 10 equally spaced

finite elements.

Although the analytical solution does not predict the

occurrence of dead core for first-order kinetics, a critical

value of Thiele modulus is determined by applying the

orthogonal collocation on finite elements. This is due to the

numerical error of the method itself. However, as shown in

Fig. 4, for orthogonal collocation on finite elements, the

calculated values of substrate concentration are approxi-

mately zero inside the predicted dead core and are within

the numerical error.

Figure 5 shows the behavior of effectiveness factor as a

function of the Thiele modulus for classical geometries

(slb.-slab, cyl.-cylinder, sph.-sphere) and reaction kinetics

of zero and first orders, calculated from the method of

orthogonal collocation on finite elements.

According to the obtained results, the orthogonal col-

location on finite elements was able to attain reliable values

for both the internal effectiveness factor and substrate

concentration profiles. This result indicates that the

orthogonal collocation method on finite elements is able to

solve similar BVPs to calculate the internal effectiveness

factor for nonlinear kinetic expressions, such as the

Michaelis–Menten kinetics and other derived expressions

incorporating inhibition terms.

Fig. 5 Effectiveness factor versus Thiele modulus for: analytical

(an.) and numerical (num.) methods; slab (slb.), cylinder (cyl.), and

sphere (sph.) geometries; zero-order (0) and first-order (1) kinetics

Table 4 Characteristics and capabilities of the numerical methods applied to solve the proposed BVPs

Shooting methods Collocation methods

Adams Stiff RK RK-4th bvode function Orthogonal Orthogonal on

finite elements

Zero-order kinetics (for low /) Solves Solves Solves Solves Solves Solves Solves

Zero-order kinetics (for high /) Solves Solves Solves Solves Solves Solves Solves

First-order kinetics (for low /) Solves Solves Solves Solves Solves Solves Solves

First-order kinetics (for high /) Diverges Diverges Diverges Diverges Inaccurate Inaccurate Solves

Initial estimate Without Without Without Without Unnecessary Necessary Necessary

Convergence rate High Medium High High Medium Medium Low

Complexity of programming on Scilab Low Low Low Low Medium High High

1724 Bioprocess Biosyst Eng (2016) 39:1717–1727

123



Assessment of numerical methods

Table 4 presents an assessment of characteristics and

capabilities of the numerical methods used to solve the

BVPs proposed in this study.

Table 4 shows that, among the various methods

employed, the only one able to solve all of the proposed

BVPs was the orthogonal collocation on finite elements

method. However, for low values of the Thiele modulus, it

may well be preferable to use the shooting method due to

its programming simplicity and high convergence rate.

Analysis of dead-core occurrence for the Michaelis–

Menten kinetics

Table 5 presents the results for orthogonal collocation

methods on finite elements to solve the BVP concerning

Michaelis–Menten kinetics.

According to Table 5, the orthogonal collocation on

finite elements was able to determinate the parameters a,

/crit, and g for all particle geometries (slab, cylinder, and

sphere) and any values of b. The calculated values of both

a and /crit indicate that the dead-core occurrence for the

Table 5 Results for orthogonal

collocation methods on finite

elements to solve the BVP

concerning Michaelis–Menten

kinetics

b = 10-5 b = 10-3

a / /crit a g a / /crit a g

1 0.5 1.4322 0 1.0000 1 0.5 2.7515 0 0.9999

1 1 1.4322 0 1.0000 1 1 2.7515 0 0.9994

1 2 1.4322 0.2826 0.7070 1 2 2.7515 0 0.7050

1 4 1.4322 0.6413 0.3535 1 4 2.7515 0.3121 0.3525

1 8 1.4322 0.8207 0.1768 1 8 2.7515 0.6561 0.1763

2 0.5 1.0009 0 1.0000 2 0.5 1.4245 0 0.9998

2 1 1.0009 0 0.9999 2 1 1.4245 0 0.9964

2 2 1.0009 0.6133 0.6217 2 2 1.4245 0.3034 0.6159

2 4 1.0009 0.8147 0.3320 2 4 1.4245 0.6547 0.3311

2 8 1.0009 0.9089 0.1715 2 8 1.4245 0.8277 0.1710

3 0.5 0.8167 0 1.0000 3 0.5 0.9663 0 0.9998

3 1 0.8167 0.3843 0.9422 3 1 0.9663 0.0392 0.9400

3 2 0.8167 0.7372 0.5933 3 2 0.9663 0.5358 0.5918

3 4 0.8167 0.8755 0.3255 3 4 0.9663 0.7696 0.3245

3 8 0.8167 0.9395 0.1703 3 8 0.9663 0.8851 0.1693

b = 100 b = 102

a / /crit a g a / /crit a g

1 0.5 61.2777 0 0.9584 1 0.5 86.5328 0 0.9249

1 1 61.2777 0 0.8397 1 1 86.5328 0 0.7628

1 2 61.2777 0 0.5427 1 2 86.5328 0 0.4829

1 4 61.2777 0 0.2770 1 4 86.5328 0 0.2503

1 8 61.2777 0 0.1385 1 8 86.5328 0 0.1252

2 0.5 31.0911 0 0.9379 2 0.5 43.7509 0 0.8935

2 1 31.0911 0 0.7743 2 1 43.7509 0 0.6989

2 2 31.0911 0 0. 4812 2 2 43.7509 0 0.4325

2 4 31.0911 0 0.2593 2 4 43.7509 0 0.2342

2 8 31.0911 0 0.1342 2 8 43.7509 0 0.1212

3 0.5 20.9338 0 0.9259 3 0.5 29.4307 0 0.8771

3 1 20.9338 0 0.7439 3 1 29.4307 0 0.6727

3 2 20.9338 0 0.4637 3 2 29.4307 0 0.4174

3 4 20.9338 0 0.2545 3 4 29.4307 0 0.2296

3 8 20.9338 0 0.1328 3 8 29.4307 0 0.1200

Bioprocess Biosyst Eng (2016) 39:1717–1727 1725

123



Michaellis–Menten kinetics is possible when verified cer-

tain conditions of diffusion–reaction within the catalytic

particle.

It was also found that for values of b close to zero, the

parameters a, /crit, and g were close to those obtained for

zero-order kinetics. For intermediate values of b, these

parameters have values between those obtained for zero-

and first-order kinetics, while for high values of b, the

calculated parameters are close to those obtained for first-

order kinetics. These results were expected since the value

of the constant K, which determines the behavior of

Michaelis–Menten kinetics between zero- and first-order

kinetics, is given by K = bs0.
The simulated profiles of substrate concentration (s) in-

side the particle for the Michaelis–Menten kinetics, a = 3

(spherical geometry) and / = 1 are shown in Fig. 6.

Figure 6 shows that the simulated profiles of s are

between the profiles for first- and zero-order kinetics.

These result agree with the asymptotic analysis of Eq. 10

presented in the mathematical modeling section. For

diffusion–reaction conditions simulated in Fig. 6, it can be

observed that the dead core occurs for low values of b (e.g.,

b B 0.01).

Effectiveness factor for the Michaelis–Menten

kinetics

Using the orthogonal collocation method on finite ele-

ments, it was possible to calculate the effectiveness factor

for three geometries and for any values of the Thiele

modulus and dimensionless Michaelis–Menten constant

(b), as shown in Fig. 7.

Figure 7 shows that, for high values of b, the calculated

values of the effectiveness factor are close to those calcu-

lated for first-order kinetics, whereas for low b values, they

are close to those calculated for zero-order kinetics. This

result was expected, given that the Michaelis–Menten

kinetics incorporates the behavior of those kinetic equa-

tions. These results confirm the validity of the orthogonal

collocation method on finite elements to solve the bound-

ary value problem associated with Michaelis–Menten

Fig. 6 Simulated profiles of substrate concentration (s) inside the

particle for the Michaelis–Menten kinetics, a = 3 (spherical geom-

etry) and / = 1

Fig. 7 Effectiveness factor versus Thiele modulus for the Michaelis–

Menten kinetics

1726 Bioprocess Biosyst Eng (2016) 39:1717–1727

123



kinetics, since that the numerical solutions were identical

to the analytical solutions for the limit cases.

Conclusions

According to the present study’s results, the following

conclusions can be drawn:

• For low values of the Thiele modulus, the use of

shooting methods is preferred due to its programming

simplicity and high convergence rate;

• For high values of the Thiele modulus, the shooting

methods can be divergent;

• For high values of the Thiele modulus, the orthogonal

collocation method and that implemented by the bvode

function of Scilab can result in inaccurate values of the

effectiveness factor;

• Among the methods studied in this article, the orthog-

onal collocation method on finite elements was the

most effective and promising to solve nonlinear

boundary value problems associated with the modeling

of diffusion–reaction processes in catalytic particles

containing immobilized enzymes;

• The orthogonal collocation method on finite elements

was successfully applied to the case of the Michaelis–

Menten kinetics, suggesting that this method can be

extended to other types of kinetic expressions, such as

those that incorporate inhibition effects;

• For the Michaelis–Menten kinetics, the dead core can

occur when verified certain conditions of diffusion–

reaction within the catalytic particle;

• The application of the concepts and methods presented

in this study will allow for a more generalized analysis

and more accurate designs of heterogeneous enzymatic

reactors.

References

1. Šekuljica NŽ, Prlainović NŽ, Jovanović JR, Stefanović AB,

Djokić VR, Dušan ŽM, Knežević-Jugović ZD (2016) Immobi-

lization of horseradish peroxidase onto kaolin. Bioproc Biosyst

Eng. doi:10.1007/s00449-015-1529-x

2. Edet E, Ntekpe M, Omereji S (2013) Current trend in enzyme

immobilization: a review. Int J Mod Biochem 2:31–49

3. Illanes A (2008) Enzyme biocatalysis: principles and applica-

tions. Springer, Netherlands

4. Bailey JE, Ollis DF (1986) Biochemical engineering fundamen-

tals. McGraw-Hill, New York

5. Doran PM (2013) Bioprocess engineering principles. Academic

Press, Waltham

6. Fogler HS (2005) Elements of chemical reaction engineering.

Prentice Hall PTR, Upper Saddle River

7. Oliveira SC (1999) Evaluation of effectiveness factor for

immobilized enzymes using Runge–Kutta–Gill method: how to

solve mathematical undetermination at particle center point?

Bioproc Eng. doi:10.1007/s004490050579

8. Granato MA, Queiroz LC (2003) Dead core in porous catalysts:

modeling and simulation of a case problem using Mathematica.

In: Proceedings of the workshop modelling and simulation in

chemical engineering, CIM–Centro Internacional de Matemática,

Coimbra, 30 June–4 July 2003

9. Bird RB, Stewart WE, Lightfoot EN (2002) Transport phenom-

ena. Wiley, New York

10. Froment GF, Bischoff KB, De Wilde J (2011) Chemical reactor

analysis and design. Wiley, Hoboken

11. Mahalakshmi M, Hariharan G (2016) An efficient Chebyshev

wavelet based analytical algorithm to steady state reaction–dif-

fusion models arising in mathematical chemistry. J Math Chem.

doi:10.1007/s10910-015-0560-0

12. Devi MR, Sevukaperumal S, Rajendran L (2015) Non-linear

reaction diffusion equation with Michaelis–Menten kinetics and

Adomian decomposition method. Appl Math. doi:10.5923/j.am.

20150501.04

13. Shanthi D, Ananthaswamy V, Rajendran L (2013) Analysis of

non-linear reaction-diffusion processes with Michaelis–Menten

kinetics by a new Homotopy perturbation method. Nat Sci.

doi:10.4236/ns.2013.59128

14. Joy RA, Meena A, Loghambal S, Rajendran L (2011) A two-

parameter mathematical model for immobilized enzymes and

Homotopy analysis method. Nat Sci. doi:10.4236/ns.2011.37078

15. Hindmarsh AC (1983) In: Stepleman RS et al (eds) Scientific

computing. IMACS Transactions on Scientific Computation,

Amsterdam

16. Brown PN, Hindmarsh AC (1989) Reduced storage matrix

methods in stiff ODE systems. Appl Math Comput. doi:10.1016/

0096-3003(89)90110-0

17. Description and use of LSODE, the Livermore solver for ordinary

differential equations (1993) NASA, Livermore. https://computa

tion.llnl.gov/casc/nsde/pubs/u113855.pdf. Accessed 16 Feb 2016

18. Shampine LF (1996) Some practical Runge–Kutta formulas.

Math Comput. doi:10.2307/2008219

19. Shampine LF, Watts HA (1980) DEPAC–design of a user ori-

ented package of ODE solvers. Sandia National Laboratories,

Albuquerque

20. Fehlberg E (1970) Klassische Runge–Kutta-Formeln vierter und

niedrigerer Ordnung mit Schrittweiten-Kontrolle und ihre

Anwendung auf Wärmeleitungs probleme. Computing. doi:10.

1007/BF02241732

21. Ascher U, Christiansen J, Russell RD (1981) Collocation soft-

ware for boundary-value ODEs. ACM T Math Softw. doi:10.

1145/355945.355950

22. Bader G, Ascher U (1987) A new basis implementation for a

mixed order boundary value ode solver. SIAM J Sci Stat Comp.

doi:10.1137/0908047

23. Ascher U, Christiansen J, Russell RD (1979) A collocation solver
for mixed order systems of boundary value problems. Math

Comput. doi:10.1090/S0025-5718-1979-0521281-7

24. Deboor C, Weiss R (1980) SOLVEBLOK: a package for solving

almost block diagonal linear systems. ACM T Math Softw.

doi:10.1145/355873.355880

25. Villadsen J, Michelsen ML (1978) Solution of differential

equation models by polynomial approximation. Prentice-Hall,

Englewood Cliffs

26. Wylie CR, Barret LC (1982) In: Wylie CR, Barret LC (eds)

Advanced engineering mathematics. McGraw-Hill, Singapore

Bioprocess Biosyst Eng (2016) 39:1717–1727 1727

123

http://dx.doi.org/10.1007/s00449-015-1529-x
http://dx.doi.org/10.1007/s004490050579
http://dx.doi.org/10.1007/s10910-015-0560-0
http://dx.doi.org/10.5923/j.am.20150501.04
http://dx.doi.org/10.5923/j.am.20150501.04
http://dx.doi.org/10.4236/ns.2013.59128
http://dx.doi.org/10.4236/ns.2011.37078
http://dx.doi.org/10.1016/0096-3003(89)90110-0
http://dx.doi.org/10.1016/0096-3003(89)90110-0
https://computation.llnl.gov/casc/nsde/pubs/u113855.pdf
https://computation.llnl.gov/casc/nsde/pubs/u113855.pdf
http://dx.doi.org/10.2307/2008219
http://dx.doi.org/10.1007/BF02241732
http://dx.doi.org/10.1007/BF02241732
http://dx.doi.org/10.1145/355945.355950
http://dx.doi.org/10.1145/355945.355950
http://dx.doi.org/10.1137/0908047
http://dx.doi.org/10.1090/S0025-5718-1979-0521281-7
http://dx.doi.org/10.1145/355873.355880

	Occurrence of dead core in catalytic particles containing immobilized enzymes: analysis for the Michaelis--Menten kinetics and assessment of numerical methods
	Abstract
	Introduction
	Mathematical modeling
	Numerical methodology

	Results and discussion
	Shooting methods
	Collocation methods
	Assessment of numerical methods
	Analysis of dead-core occurrence for the Michaelis--Menten kinetics
	Effectiveness factor for the Michaelis--Menten kinetics

	Conclusions
	References