CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 70, 2018 

A publication of 

 

The Italian Association 
of Chemical Engineering 
Online at www.aidic.it/cet 

Guest Editors: Timothy G. Walmsley, Petar S. Varbanov, Rongxin Su, Jiří J. Klemeš 
Copyright © 2018, AIDIC Servizi S.r.l. 

ISBN 978-88-95608-67-9; ISSN 2283-9216 

Research and Mathematical Modelling of Direct 

Bioelectrocatalytic Oxygen Reduction by Laccase 

Violetta A. Vasilenkoa*, Irina N. Arkadevaa, Vera A. Bogdanovskayab, Evgenia A. 

Fokinaa, Eleonora M. Koltsovaa 

aD.Mendeleev University of Chemical Technology of Russia, Miusskaya sq., 9, Moscow, 125047, Russia  
bA.N. Frumkin Institute of Physical chemistry and Electrochemistry Russian Academy of Sciences, Leninsky prospect, 31, 

 Moscow, 199071 Russia 

 kolts@muctr.ru 

At this present work, the process of direct bioelectrocatalytic oxygen reduction by laccase was investigated 

using electrochemical methods and mathematical modelling. The maximum of achieved current density is 640 

μA cm-2. The developed mathematical model includes the mass and potential conservation equations and takes 

into account the porous structure of the electrode active layer. The equations of the mathematical model were 

solved by numerical methods with own developed software package. The profiles of the components 

concentration on time and thickness of the active layer of the electrode were obtained. The adequacy of the 

mathematical model was tested by comparing the experimental and calculated values of the electrodes activity 

at different loadings of the carbon material. By the golden ratio method, the optimal electrode carbon material 

loading of 0.92 mg cm-2 was found. 

1. Introduction 

Presently the primary interest of researchers in the field of mathematical modelling of fuel cells (FCs) is related 

to processes of traditional fuel cells (low-temperature with proton-conductive polymer electrolyte and high-

temperature solid oxide), while the papers with mathematical modelling of biofuel cells (BFC) are much less. 

This area of research has become especially accessible only in the last decade, and the published to date works 

on this topic are general, seeking to cover all operational aspects of the BFC.  

Mathematical modelling and simulation can play a fundamental role in the further understanding and 

development of BFC. Adequate mathematical models cooperated with laboratory studies, can be used for 

studying the reaction medium, engineering of new electrode designs and accelerate the development of 

practical-useful systems. According to the review published by Rajendran et al. (2017) most mathematical 

models of biofuel cells are based on the reaction-diffusion differential equations. The primary type of models is 

based on the principles of formal kinetics. The major part of the models predicts the characteristics of a BFC 

with separation membrane and with the usage of electron transfer mediators. But practically there is no a 

mathematical model describing the constructions of membraneless and mediatorless BFC under the direct 

electron transfer of current generating reactions.  

There are some investigations where mathematical modelling and calculations were carried out with the help of 

commercial package software, e.g. CFD-ACE (Bedekar et al., 2007), COMSOL Multiphysics (Chan et al., 2012), 

ANSYS FLUENT (Mitrichev et al., 2016), FEM-LAB (Kjeang et al., 2006) or use for calculations the particular 

software (MATLAB) (Osman et al., 2013). These models operate with the enzyme surface concentration and 

take it equal to the total amount of enzyme adsorbed at the electrode. Using these models, it is impossible to 

determine the factors which can affect the final characteristics of the electrodes. In the previous study (Arkadeva 

et al., 2017) we calculated the laccase distribution within the active layer of the electrode as a result of 

spontaneous adsorption immobilization. These data have been taken in the present work as initial conditions 

for determining the electrocatalytic activity of the electrode as a function of the thickness of the active layer. 

                                

 
 

 

 
   

                                                  
DOI: 10.3303/CET1870269 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Please cite this article as: Vasilenko V.A., Arkadeva I.N., Bogdanovskaya V.A., Fokina E.A., Koltsova E.M., 2018, Research and mathematical 
modelling of direct bioelectrocatalytic oxygen reduction by laccase , Chemical Engineering Transactions, 70, 1609-1614  
DOI:10.3303/CET1870269   

1609



2. Experimental and model 

2.1 Chemicals and methods 

In this study, the Laccase Trametes Versicolor - copper-containing enzyme, which can directly reduce (without 

the use of any mediator) oxygen to water was used. It was produced by the procedure described by Gorshina 

et al. (2006). Experimental researches were carried out according to Bogdanovskaya et al. (2017). 

Electrochemical measurements were conducted in a three-electrode cell with separated electrode spaces. The 

floating electrode was used as a cathode, and it was made in the form of the tablet from hydrophobized carbon 

black, electrode thickness was equal to 2 mm, the surface area of contact with electrolyte was equal to 1 cm2. 

Some amount of multiwalled carbon nanotubes provided by D.I. Mendeleev University of Chemical Technology 

of Russia (CNT, specific surface area measured by BET method SBET=210 m2 g-1) was pressed onto the 

electrode surface. The procedure of CNT fabrication was described by Bogdanovskaya et al. (2016). Then the 

electrode was placed on the surface of laccase solution, which contained 0,076 mg ml-1 of laccase. Adsorptive 

spontaneous immobilization was carried out for 2 hours. The quantity of adsorbed enzyme was determined 

spectrophotometrically using a preconstruction calibration curve as the difference of enzyme content in the 

solution before and after adsorption process. 

A platinum wire was used as a counter electrode, and a saturated Ag/AgCl electrode was used as a reference 

electrode. Studies were carried out in 0,2 M phosphate-acetate buffer solution (pH = 4,5). Polarization curves 

were registered in the oxygen atmosphere using the potentiostat IPC-Pro 3A. 

2.2 Mathematical simulation 

The developed mathematical model assumes two modelling domains (Figure 1). The first one is the volume of 

the electrolyte (0<x<δ) and the second one is the active layer, which is the porous material filled with electrolyte 

(δ<x<h). 

The area of the electrolyte volume can be described by the ideal mixing model, i.e., there are no concentration 

gradients in this region. Since the process proceeds in a buffer solution with a constant supply of oxygen, the 

concentration values of the protons and dissolved oxygen in this region remain constant. Furthermore, in this 

region does not occur electrochemical reaction. 

 

 

Figure 1: Schematic representation of modelling domains 

In the active layer there is the oxygen bioelectroreduction by laccase (L), which schematically can be 

represented as: 

22

1

1

LOOL
k

k

  (1) 

ОHLНеLO
ck

22 244    (2) 

It includes the reversible stage of enzyme-substrate complex LO2 production – reaction (1) and the electron 

transfer with dissociation of the enzyme-substrate complex to product – reaction (2). 

The rates of these electrode reactions can’t be described by the Michaelis-Menten equation because only the 

first stage is enzymatic. The second stage is electrochemical and can described by Buttler-Volmer equation. 

The rates of the reactions (1) and (2) respectively (the subscript «с» characterizes the cathodic process): 

211 OLcckW   (3) 

211 LOckW    (4) 

1610













RT

nF
expckW LOcc


 2

 (5) 

where cj, mol m-3 – concentration of the components; α – charge transfer coefficient; n – number of electrons in 

the electrochemical reaction; ηс , V – the potential drop between the carbon material (CM) of the active layer of 

the electrode and the electrolyte; k1γ, k-1γ, kcγ – modified reaction rate constants, related to the corresponding 

rate constants by: k1γ=k1
γkγ, m3 mol-1 s-γ; k-1γ=k-1

γ, s-γ; kcγ=kc
γ, s-γ. The γ subscript characterizes the process, 

which takes place within the porous CM layer, and the γ superscript corresponds the material’s porosity 

(determines the part of pores, which are opened for electrolyte flow), 0<γ<1 (Koltsova et al., 2000); F, С mol-1 – 

Faraday constant; R, J mol-1 K-1 – universal gas constant; Т, К – temperature. 

It was assumed that the enzyme molecules are presented only within the active layer of the electrode, and they 

are firmly adsorbed on the surface of the CM and can't pass into solution. The amount of the enzyme and 

enzyme-substrate complex is related to the total amount of the laccase adsorbed on the electrode by the 

following equation: 

)x,t(LO)x,t(LL
ccc

2
0   (6) 

where 𝑐𝐿0, mol m-3 – the initial concentration of laccase, distributed within the active layer (obtained in Arkadeva 

et al., 2017). 

The flow of dissolved oxygen and protons in the electrolyte arises from their migration in the electric field (for 

protons) and diffusion due to the concentration gradient, and it can be described by the Nernst-Planck equation. 

Since the active layer is CM, which, as expected, due to its roughness, porosity and other irregularities can be 

described using a mathematical apparatus developed for fractal structures. Thus, the equation of material 

balance for an element of a small volume of the active layer, as in Koltsova et al. (2000), has the form: 

j
j

j
lj

jj
j

S
x

с
D

x

Ф

x

с

RT

F
Dz

t

с




















2

2





 (7) 

where zj – particle’s charge; Фl , V – the potential of the liquid (ion-conducting) phase; Sj – j-particles source 

term (the multiplication of the stoichiometric coefficient of the particle by the reaction rate with the corresponding 

sign – «-» for reagents, «+» for products); t, s –time; Djγ, m2 s-γ– the component’s diffusion coefficient within the 

porous active layer, related to general diffusion coefficient by Djγ=Dj
γS(1-γ); S, m – channel cross-sectional area 

available for the electrolyte flow; х, m – electrode thickness coordinate; Dj, m2 s-1 – particle’s diffusion coefficient. 

According to the Eq(7), the equations for the components concentration change within the active layer of the 

cathode are: 

112

2
222










WW

x

c
D

t

c OОO




 (8) 

c
H

H
lHH

H
H W

x

c
D

x

Ф

x

c

RT

F
Dz

t

c
4

2

2


















 













 (9) 

c
L WWW

t

c





11



 (10) 

The potential drop was determined as: 

η = Фs – Фl - Voc (11) 

where the Voc – the open circuit voltage. 

The connection between the potential and the current density was expressed by the Poisson equation and for 

the potential of the solid phase Фs (CM, electron-conducting) and of the liquid phase Фl as in Do et al. (2014) 

taking into account the Bruggeman correlation for the difference in the conductivity of the ion and electronic 

phases: 

i
x

Ф
)(

x

S
S

. 
















 511  (12) 

For the potential of the ionic (liquid) phase: 

i
x

Ф

x

l
l

. 















 51  (13) 

1611



where ε – volume fraction of the liquid phase in the active layer of the electrode; κs, κl, А (V·m)-1 – the conductivity 

of the electron (s) and the ion (l) conducting phases; i, А m-3– local current density. 

The local current density is related to the rate of the electrode reaction by the Faraday equation: 

cnFWi   (14) 

The value of the total current density was obtained by spatial integration of the local current densities along the 

thickness of the active layer of the electrode: 




0

dx)x,t(iI t  (15) 

where δ, m – active layer thickness. 

Eq(8) and Eq(9) are the differential equations of parabolic type, solved by the sweep method. Eq(10) is an 

ordinary differential equation of the first order, solved by the implicit Euler method. Eq(12) and Eq(13) are 

equations of elliptic type, solved using the establishment method, adding the time derivative. 

In contrast to works, where model equations were solved using the commercial software, in this work the model 

equations were solved by the finite-difference method, the program module was implemented in the C++ 

programming language. 

3. Results and discussion 

3.1 Results of electrochemical experiments 

Figure 2 shows the polarization curves for the oxygen electroreduction by laccase for electrodes with different 

amounts of CM applied to the electrode. It was found the limiting current of the electroreduction of oxygen by 

laccase increases with the increasing of the amount of CM, and the greatest activity of the electrode was 

observed with 0.94 mg of CM.  

 

Figure 2: Polarization curves for the oxygen electroreduction by the electrode with laccase, adsorbed on CM 

CNT, which was pressed into the electrode with an amount: 1 – 0,37 mg; 2 – 0,46 mg; 3 – 0,48 mg; 4 – 0,73 

mg; 5 - 0,86 mg; 6 – 0,94 mg. Potential sweep rate 1 mV s-1 

In Arkadeva et al., 2017 it was shown that further increase of the CM amount leads to decreasing the activity of 

the laccase-based electrode. The efficiency of the enzyme adsorption also decreases. In Bogdanovskaya et al. 

(2017) it was assumed that usage of this type of CM allows laccase to adsorb not only on the surface of the 

material, but penetrate into its mesopores and also other voids. If the optimal CM loading for a given electrode 

is overloaded, the proportion of voids decreases, and accordingly the size of the potentially suitable space for 

laccase adsorption. 

The obtained data were applyed to develop a mathematical model of the bioelectrocatalytic oxygen reduction 

by laccase. 

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0

350400450500550600650700750

C
u
rr

e
n
t 

d
e
n
s
it
y
, 
m

A
 c

m
-2

Potential, mV

6
5

4
3

2
1

1612



3.2 Results of mathematical modeling 

The solution of the mathematical equations was carried out for experiments with the loading of CM 0,6 mg. The 

reaction constants and the conductivity coefficients were found by the scanning method. The simulation 

parameters are listed in the Table 1. The literature data were taken as the initial approximations. 

Table 1: Simulation parameters 

Parameter  Value Reference Parameter Value Reference 

T, K 298 experimental kc, m3 mol-1 s-γ 1,43·10-13 fitted 

cL
0

(t=0, x) mol m-3 1,375·10-3 Arkadeva et al. (2017) κl, A m-1 V-1 0,055 fitted, initial Chan et. 

al. (2012) 

cO2
bulk, mol m-3 1,38 experimental κs, A m-1 V-1 0,046 fitted, initial Chan et. 

al. (2012) 

cH+
bulk, mol m-3 31,6·10-3 experimental ε 0,8 fitted, initial Osman 

et al. (2013, b)  

DO2, m2 s-1 2·10-9 Osman et al. (2013, a) α 0,5 experimental 

DH+, m2 s-1 1·10-9 Chan et al. (2012) n 4  

k1, m3 mol-1 s-γ 1000 fitted, initial Do et al. 

(2014) 

Voc, V 0,18 experimental 

k-1, s-γ 646 fitted, initial Do et al. 

(2014) 

 0,94  Arkadeva et al. 

(2017) 

 

The equations of the mathematical model were calculated using the parameters from the Table 1 and the values 

of the current density as a function of the cell voltage were obtained. The model was tested for different CM 

loading: the Figure 3 shows the experimental and calculated data for the values of current density depending 

on the mass of the CM. The results of the calculation are consistent with the experimental data. By the usage 

of golden ratio method it was found that the optimal electrode carbon material loading is 0.92 mg cm-2. 

 

Figure 3: The experimental and the calculated volumes of the electrode activity on the loading of CM 

The calculated distributions of the components over the thickness of the electrode active layer in 100 s after the 

start of the bioelectrocatalytic oxygen reduction are shown on Figure 4. 

  
a b 

Figure 4: Calculated concentration distributions in the active layer: a) laccase and enzyme-substrate complex; 

b) oxygen and protons  

-0.7

-0.65

-0.6

-0.55

-0.5

0.55 0.65 0.75 0.85 0.95 1.05 1.15

C
u
rr

e
n
t 

d
e
n
s
it
y
, 
m

A
 c

m
-2

Mass of CM, mg

Experimental Numerical

0.4

0.6

0.8

1.0

0 1 2 3 4 5 6 7

C
o
n
c
e
n
tr

a
ti
o

n
, 

m
m

o
l m

-3

CM layer thickness, μm
LO2 L

0

10

20

30

40

0.8

0.9

1.0

1.1

1.2

1.3

1.4

0 1 2 3 4 5 6 7 P
ro

to
n
 c

o
n
c
e
n
tr

a
ti
o

n
,

m
m

o
l m

-3

O
x
y
g
e
n
 c

o
n
c
e
n
tr

a
ti
o

n
, 

m
o
l 
m

-3

CM layer thickness, μm

O2 H+

1613



It was revealed that at the first steps of the calculation the components concentrations varied significantly, and 

then remained nearly unchanged. Probably such character of the dependence is due to the high values of the 

reactions rate constants.  

In the future, the obtained data will be used to optimize the design of the enzymatic FC cathode. 

4. Conclusions 

In this paper, a detailed dynamic mathematical model was developed on the basis of the experimental studies 

of the activity of the electrode based on laccase in the bioelectroreduction of oxygen to water. As a result of the 

simulation, the distribution profiles of the concentrations of the components along the thickness of the active 

layer of the cathode were obtained. The dependences of the current density on the cell voltage for different 

loads of the carbon material (CNT) were calculated. The calculated data are consistent with the experimental 

data, which indicates the adequacy of the model. The optimal loading of the electrode by the carbon material 

was also found. 

Acknowledgments 

Authors gratefully acknowledge Doctor of chemical sciences, Prof. Gorshina N.N. and Doctor of chemical 

sciences, Prof. Yaropolov A.I. for providing the laccase and research assistance. 

This study is partially supported by the state contract № 13.10932.2018/10.11. 

References 

Arkadeva I.N., Bogdanovskaya V.A., Vasilenko V.A., Fokina E.A., Koltsova E.M., 2017, Laccase spontaneous 

adsorption immobilization: experimental studies and mathematical modeling at enzymatic fuel cell cathode 

construction, IOP Conf. Series: Earth and Environmental Science, 83, 127-134. 

Bedekar A.S., Feng J.J., Krishnamoorthy S., Lim K.G., Palmore G.T.R., Sundaram S., 2007, Oxygen limitation 

in microfluidic biofuel cells, Chemical Engineering Communication., 195, 256-266. 

Bogdanovskaya V.A., Arkadeva I.N., Osina M.A., 2017, Bioelectrocatalytic oxygen reduction by laccase 

immobilized on various carbon carriers, Russian Journal of Electrochemistry, 53 (12), 1352-1362. 

Bogdanovskaya V. A., Koltsova E. M., Zhutaeva G. V., Radina M. V., Kazanskii L. P., Tarasevicha M. R., 

Skichko E. A., Gavrilova N. N., 2016, Physico- chemical properties of carbon nanotubes as supports for 

cathode catalysts of fuel cells. Surface structure and corrosion resistance, Protection of Metals and Physical 

Chemistry of Surfaces, 52, 45-54. 

Chan D.S., Dai D.J., Wu H.S., 2012, Dynamic modeling of anode function in enzyme-based biofuel cells using 

high mediator concentration, Energies, 5, 2524-2544. 

Do T.Q.N., Varničić M., Hanke-Rauschenbach R., Vidaković-Koch T., Sundmacher K., 2014, Mathematical 

modeling of a porous enzymatic electrode with direct electron transfer mechanism, Electrochimica Acta, 

137, 616-626. 

Gorshina E.S., Rusinova T.V., Biryukov V.V., Morozova O.V., Shleev S.V., Yaropolov A.I., 2006, The dynamics 

of oxidase activity during cultivation of basidiomycetes from the genus Trametes Fr., Applyed Biochemistry 

and Microbioljgy, 42, 638-644.  

Kjeang E., Sinton D., Harrington D.A., 2006, Strategic enzyme patterning for microfluidic biofuel cells, Journal 

of Power Sources, 158, 1–12. 

Koltsova E.M., Vasilenko V.A., Tarasov V.V., 2000, Numerical methods for solving transport equations in fractal 

media, Russian Journal of Physical Chemistry A, 74 (5), 848-850. 

Mitrichev I., Vasilenko V., Scherbakov A., Koltsova E., 2016, Simulation of PEMFC with Pt/CNT nanostructured 

cathode catalyst, Chemical Engineering Transactions, 52, 1003-1008. 

Osman M. H., Shah A. A., Wills R. G. A., Walsh F. C., 2013a, Mathematical modelling of an enzymatic fuel cell 

with an air-breath cathode, Electrochim. Acta, 112, 386-393. 

Osman M. H., Shah A.A., Wills R. G. A, 2013b, Detailed mathematical model of an enzymatic fuel cell, Journal 

of The Electrochemical Society, 160 (8), F806-F814. 

Rajendran L., Kirthiga M., Laborda E., 2017, Mathematical modeling of nonlinear reaction–diffusion processes 

in enzymatic biofuel cells, Current Opinion in Electrochemistry, 1, 121 – 132. 

Song Y., Penmatsa V., Wang C., 2014, Modeling and simulation of enzymatic biofuel cells with three-

dimensional microelectrodes, Energies, 7, 4694-4709. 

 

1614

https://www.sciencedirect.com/science/article/pii/S2451910316300187#!
https://www.sciencedirect.com/science/article/pii/S2451910316300187#!
https://www.sciencedirect.com/science/article/pii/S2451910316300187#!