*Article* **Continuum Scale Non Newtonian Particle Transport Model for Hæmorheology**

**Torsten Schenkel 1,\* and Ian Halliday 2**

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** We present a continuum scale particle transport model for red blood cells following collision arguments, in a diffusive flux formulation. The model is implemented in FOAM, in a framework suitable for haemodynamics simulations and adapted to multi-scaling. Specifically, the framework we present is able to ingest transport coefficient models to be derived, prospectively, from complimentary but independent meso-scale simulations. For present purposes, we consider modern semi-mechanistic rheology models, which we implement and test as proxies for such data. The model is verified against a known analytical solution and shows excellent agreement for high quality meshes and good agreement for typical meshes as used in vascular flow simulations. Simulation results for different size and time scales show that migration of red blood cells does occur on physiologically relevany timescales on small vessels below 1 mm and that the haematocrit concentration modulates the non-Newtonian viscosity. This model forms part of a multi-scale approach to haemorheology and model parameters will be derived from meso-scale simulations using multi-component Lattice Boltzmann methods. The code, haemoFoam, is made available for interested researchers.

**Citation:** Schenkel, T.; Halliday, I. Non Newtonian Particle Transport Model for Hæmorheology. *Mathematics* **2021**, *9*, 2100. https:// doi.org/10.3390/math9172100

Academic Editors: Victor Mitrana and Alexander Zeifman

Received: 29 May 2021 Accepted: 25 August 2021 Published: 30 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Keywords:** haemorheology; blood flow modelling; particle transport; numerical fluid mechanics

#### **1. Introduction**

Blood is a non-Newtonian fluid with very complex behaviour deriving from a mesoscopic composition which—minimally described—is a dense, mono-disperse suspension of deformable vesicles suspended in incompressible plasma. Accordingly, blood rheology is dominated by the interaction of cells, with a multitude of models having been proposed to account for such meso-scale effects as deformation, aggregation, and rouleaux formation which underline emergent macroscopic flow properties like concentration dependant viscosity and shear thinning. The authors are currently developing a multi-scale approach, explicitly modelling meso-scale effects using Lattice Boltzmann Models (LBM), in which erythrocyte mechanics are fully resolved, while describing the macro-scale rheology using particle transport modelling and quasi-mechanistic non-Newtonian rheology models described using traditional Eulerian numerics. We give the rationale of this pairing of methodologies shortly, in the next sub-section. In essence, though, it is intended that the latter model will eventually be parameterised using LBM data. Here, we present the continuum mechanical part of this multi-scale modelling approach, which is designed to facilitate the simulation of realistic vessel geometries and spatially complex flow patterns.

Multi-scale models are necessary in many applications. In meteorology, the Eulerian grids of the macroscopic simulation are characterised by a 1 km spatial resolution. Clearly it is necessary to account for more local variation of field variables and macroscopic atmospheric simulations are parameterised based on data gathered from, e.g., semi-analytic meso-scale models of convection [1]. (Note, we use the term meso-scale (macroscale) to identify the shortest (largest) length scale considered. We use the term continuum-scale to

describe a scale where particle dynamics are neglected and integral bulk material parameters are used—within the scope of this paper macroscale is treated as continuum-scale).

Similarly, in the simulation of turbulent flow fields, temporal or spatial averages are used to derive macroscopic turbulence properties, like an isotropic turbulent viscosity, which describes the statistical effect of the inherently chaotic, and anisotropic nature of meso-scale and micro-scale turbulence. While the micro-scale can be simulated (Direct Numerical Simulation), this approach is limited by size, Reynolds number and computational resource and thus it is mostly used in fundamental research and in the development and parametrisation the macro-scale model.

In a high particle load suspension like blood, many types of mechanical interactions between particles and carrier fluid occur [2]. A variety of models are used from 1 to d reduced order models, through meso-scale models to continuum-scale models [3]. Meso-scale modelling, e.g., Particle Dynamics Methods [4], or the multi-component Lattice Boltzmann Method [5] has a widely acknowledged facility for Lagrangian particulate flows [6–10] and has been employed to describe these interactions and the dynamics of the collisions in detail. Very recently, we have developed a single framework, three-dimensional methodology for capturing detailed, particle-scale interactions between neutrally buoyant suspended vesicles (i.e., erythrocytes), using our novel chromodynamics multi-component Lattice Boltzmann method variant [11]. We have previously shown that this same essential model is able to capture detailed hydrodynamic interactions, lubrication effects and ballistic collisions between transported particles, in two dimensional simulations containing O(1000) liquid drops at high volume faction [12]. Current work aims to scale these models up, to allow extraction, through statistical averaging and measuring mechanical dissipation, of relevant transport coefficients, e.g.,

$$\nu = \nu(\phi, \dot{\gamma}), \quad D = D(\phi, \dot{\gamma}). \tag{1}$$

Above, *ν* and *D* denote kinematic viscosity and diffusion co-efficient and *φ* and *γ*˙ denote particle density and local shear rate.

As with direct numerical simulation in turbulence modelling, finite computational resource means that these detailed explicit particulate models are limited to small volumes containing relatively few particles in their simulation domain (an the order of magnitude of hundreds to thousands at the time of writing). To address the much greater scales of medical significance, it is, therefore, necessary to develop macro or continuum scale models, which encapsulate an integrated effect of these interactions, without explicitly resolving them. Crucially, these continuum models must be amenable to parametrisation, in the present context using meso-scale information such as that encapsulated in the data of Burgin [12].

Based on the work of Leighton and Acrivos [13], and Phillips et al. [14], we present here a macro-scale, or continuum-scale model for haematocrit transport, which, together with modern formulations of "quasi-Newtonian" (i.e., not considering viscoelastic and viscoplastic effects) non-Newtonian model, allows for the simulation of macroscopic flows, accounting for the integral effect of meso-scopic phenomena.

Currently, the models which have been proposed for continuum-scale haematocrit transport models can, roughly, be divided [15] into suspension balance models [13,16] and diffusive flux models [14].

Suspension balance models use an Euler–Euler mixture modelling approach, where the carrier fluid and the particle load are represented as separate species, with a transport Equation (typically convection-diffusion) and physical transport properties for each species; in diffusive flux models, the suspension is modelled as a single species, with the particle volume fraction being modelled as a scalar property, which influences the bulk transport coefficients.

Our macroscopic model is a particle transport model after Phillips [14]. It follows the collision arguments by Leighton and Acrivios [13] by describing the particle migration based on the gradients of shear strain, concentration and viscosity. The local concentration of haematocrit is then used to establish the local effective viscosity.

The haematocrit distribution emerges, locally, from a balance of the competing effects of diffusion-driven flux down the local concentration gradient (which, for steady, developed flow, lies in the radial direction) and shear-induced migration (which, for steady, developed, flows, would be in the negative radial direction). Clearly, the problem of flow-concentration coupling is made still more complex since the diffusion co-efficient varies with local shear rate and particle concentration.

Previous attempts to implement this class of particle transport model in a finite volume code were either limited to a single viscosity model (due to the necessity of linearising the source terms [17]) or they had to introduce damping terms (to suppress instability [18], mostly due to the fact that they attempted to implement the non-linear transport equation for haematocrit into commercial software packages (Fluent), which does not allow a fully implicit additional transport equation in user subroutines).

The framework we present in this work expands on those commonly used in the simulation of clinically relevant vascular flows in the following important aspects:

It is implemented using the *Field Operation And Manipulation* (FOAM) framework. FOAM, or OpenFOAM, is an open source library which allows easy implementation of Finite Volume Method (FVM) solvers and achieves stability without artificial damping by implementing the non-linear terms in a hybrid explicit-implicit fashion which allows stable simulations even in realistic physiological geometries at high rates.

#### **2. Materials and Methods**

#### *2.1. Particle Migration Model*

In a particle transport model, the transport of haematocrit is dominated by advection —following the bulk flow—variations in concentration are evened by diffusive processes, and the migration within the bulk is driven by a migration pressure. This migration pressure is the result of two phenomena: (1) spatial variation of collision (interaction) frequency and (2) spatial variation of viscosity [14]

A detailed treatise on the rationale behind the compression arguments can be found in Leighton and Acrvios, and Phillips [13,14], we only give a brief outline at this point.

#### 2.1.1. Spatial Variation of Collision Frequency

Particles that are moving relative to each other in neighbouring shear surfaces will experience collisions. The frequency of these collisions is proportional to the shear rate *γ*˙ , the particle concentration *φ*, and the particle collision radius *a*. In a field of constant concentration and constant shear, *γφ*˙ = const, the collisions are in equilibrium either side of the shear surface, and no net migration will occur. In the presence of gradients of shear rate or concentration, the imbalance of collisions will lead to a "migration pressure" down the gradient. This collision driven migration pressure can be described as a function of *a*∇(*γφ*˙ )). Using a proportionality factor of *Kc* and assuming a displacement proportional to the particle radius *a*, the migratory flux *Nc* due to variations in collision frequency can be expressed as (using the chain rule):

$$N\_{\mathfrak{c}} = -\mathcal{K}\_{\mathfrak{c}}a^2(\phi^2 \nabla \dot{\gamma} + \phi \dot{\gamma} \nabla \phi) \tag{2}$$

#### 2.1.2. Spatial Variation of Viscosity

The displacement of particles after a collision is moderated by viscous effects. In a constant viscosity field the displacement is isotropic and thus balanced with no net migration effects. In a viscosity gradient the displacement will be less damped in direction of the lower viscosity, leading to a net migration effect down the viscosity gradient.

The displacement velocity is proportional to the relative change in viscosity over a distance that is of order *a*: *a*(1/*μ*)∇*μ*. With the displacement frequency scaling with *γφ*˙ , and a proportionality factor of *Kμ*, the migratory flux due to viscosity gradient can be described as (flux is proportional to *φ*):

$$N\_{\mu} = -K\_{\mu}\dot{\gamma}\phi^2 \left(\frac{a^2}{\mu}\nabla\mu\right) \tag{3}$$

The scalar transport equation for haematocrit, *φ*, is then (neglecting molecular diffusion, Brownian motion), where *D*/*Dt* is the total differential:

$$\frac{D\Phi}{Dt} = \begin{array}{c} \nabla \cdot (\mathbf{N}\_c + \mathbf{N}\_\mu) \\\\ \frac{D\Phi}{Dt} = \\ + \quad a^2 \mathbf{K}\_c \nabla \cdot \left(\phi^2 \nabla \dot{\gamma}\right) \\\\ + \quad a^2 \mathbf{K}\_\mu \nabla \cdot \left(\phi^2 \nabla \dot{\gamma}\right) \\\\ + \quad a^2 \mathbf{K}\_\mu \nabla \cdot \left(\dot{\gamma}\phi^2 \frac{1}{\mu} \nabla \mu\right) . \end{array} \tag{5}$$

with *a*, particle radius, *γ*˙ , shear strain rate magnitude, *μ*, dynamic viscosity, *Kc* and *Kμ*, collision parameters.

Typically, the viscosity is *μ* = *f*(*γ*˙ , *φ*), which makes the last source term non-linear, which can, in turn, make the solution of this transport equation difficult.

Previous attempts to solve this problem analytically or implement this type of migration model in a numerical model used linearisation of this source term, which involves the derivative of *μ* in both *γ*˙ and *φ*, and thus limits the model to a specific viscosity model, for which it has been implemented [17,19]. Our current implementation deals with the non-linear viscosity source term in a way that leaves the viscosity gradient term intact and is thus agnostic to the rheology model used.

#### *2.2. Rheology Models*

It is obvious from the third RHS term in Equation (5), that the particle transport strongly depends on the rheology model it is coupled with. This model implementation aims to be independent of the rheology model. The draw-back of this approach is that errors present in the rheology model, which influence the particle transport, cannot be calibrated out with the parameters of the migration model alone, but the combined set of model parameters will need to be found for any new rheology model that is to be implemented.

Typically, only the shear thinning effects are taken into account, when modelling the non-Newtonian properties of blood in CFD. Common models are of the Carreau and Casson types. In these models, the haematocrit concentration is only used as a bulk parameter in the parametrisation, if at all. Our framework, incorporating the transport of haematocrit, allows the rheological model to take the local particle concentration into account when calculating the local, effective viscosity.

The rheology models that have been implemented and tested in this study are the concentration dependent Krieger–Dougherty model [20], the Quemada model [21–23] with modification by Das [24] (and a new parameter set, which avoids the singularity problem commonly associated with this model), an extended Krieger model, accommodating shear thinning and aggregation effects [25], a Casson model with haematocrit dependence following Merril et al. [24,26], and a modified Carreau type model, proposed by Yeleswarapu [27]. All model parameters have been fitted to the experimental data of Brooks [28] (Figure 1).

#### 2.2.1. Krieger–Dougherty Model

The traditional Krieger–Dougherty model [20] was developed to describe the rheology of high volume ratio suspensions of rigid spherical particles. Rigid, spherical particles do not exhibit shear-thinning behaviour, so the Krieger–Dougherty model is only dependent on the haematocrit concentration *φ*. It shows a singularity for *φ* = *φ*∗, where *φ*∗ is the haematocrit concentration for which the suspension does stop to behave like a fluid. For rigid spheres *φ*∗ = 0.68 [20], while for blood it can go up to *φ*∗ = 0.98, which is usually attributed to the deformability of the erythrocytes [25].

$$
\mu = \mu\_P \left( 1 - \frac{\Phi}{\Phi^\*} \right)^{-n}. \tag{6}
$$

The parameter *n* = *kφ*∗ is often set to *n* = 2, but more commonly to the high shear limit of *n* = 1.82 for *φ*∗ = 0.68 [22,29], which is also the value used in this work to allow comparison with the results from Phillips and others [14,17,19]. *μP* is the Newtonian viscosity of the liquid phase (plasma).

In this study, the Krieger–Dougherty model is not used as for modelling blood viscosity but as a reference model for verification and validation.

#### 2.2.2. Quemada Model

The Quemada model is based on "optimisation of viscous dissipation" [21]. In its original form it is formulated as a Newtonian, concentration dependent viscosity:

$$
\mu = \mu\_P \left( 1 - k \phi \right)^{-2},
\tag{7}
$$

with *k* being related to the packing concentration and (for the high shear limit) given as: *k* = 2/*φ*∗. In this form it is closely related to the Krieger–Dougherty model (Equation (6)). In its non-Newtonian form *k* is expressed as [22,23]:

$$k = \frac{k\_0 + k\_\infty\sqrt{\dot{\gamma}/\dot{\gamma}\_c}}{1 + \sqrt{\dot{\gamma}/\dot{\gamma}\_c}},\tag{8}$$

where *k*0 and *k*∞ are the intrinsic viscosities at zero and infinite shear, respectively, and *γ*˙ *c* is a critical shear rate.

The shear rate magnitude *γ*˙ is defined as

$$
\gamma := \sqrt{2D : D\_\prime} \tag{9}
$$

with *D*, the symmetric part of the velocity gradient tensor.

Different parameter fits have been proposed for *k*0, *k*∞, *γ*˙ *c*. Cokelet [26,30] proposed:

$$k\_0 \quad = \; \exp(a\_0 + a\_1\phi + a\_2\phi^2 + a\_3\phi^3) \tag{10}$$

$$k\_{\infty} \quad = \; \exp(b\_0 + b\_1 \phi + b\_2 \phi^2 + b\_3 \phi^3) \tag{11}$$

$$
\gamma\_c \quad = \; \exp(c\_0 + c\_1 \phi + c\_2 \phi^2 + c\_3 \phi^3) . \tag{12}
$$

Das [24] noted that Cokelet us parameter set causes the viscosity to be non-monotonous over haematocrit concentration for low shear, and exhibits singularities for zero shear. Das changed the parameter fit for *k*0 to

$$k\_0 = a\_0 + \frac{2}{a\_1 + \phi'} \tag{13}$$

which results in a monotonous behaviour for low shear (the lowest shear measured in the Brooks dataset is around *γ*˙ = 0.15 s−1), but still shows a singularity for *φ* = 80.4%. While this is outside the haematocrit values typically encountered in clinical practice, it can still pose a problem if cell migration is taken into account, which will concentrate cells in the core region. In order to overcome this problem, a new parameter set, based on Das's formulation, is derived in this work, which does not show a singularity. Figure 2 shows viscosity over shear rate for low shear rate (*γ*˙ = 0.15 s−1) and zero shear rate. While all the curves show a good fit with the data, the new parameter set does show monotonous behaviour throughout and no singularity below the critical haematocrit.

#### 2.2.3. Modified 5 Parameter Krieger Model

Hund et al. [25] proposed and developed a quasi-mechanistic extension to the Krieger– Dougherty model.

Starting from the traditional formulation of the Krieger–Dougherty model:

$$
\mu = \mu\_P \left( 1 - \frac{\phi}{\phi^\*} \right)^{-n},
\tag{14}
$$

describing the haematocrit dependence, the shear-thinning behaviour is introduced by a variable exponent *n*:

$$n = n\_{\infty} + \begin{cases} 0, \,\phi < \phi\_{st} \\ n\_{st}, \phi > \phi\_{st} \end{cases} \tag{15}$$

where *φst* is the threshold haematocrit concentration below which no shear-thinning is observed. Based on Brooks [28], this threshold is around *φ* = 0.15, and *n*∞ is modelled using a exponential dependency on *φ*:

$$m\_{\infty} = a + b \exp(-c \,\phi). \tag{16}$$

Hund's [25] shear-thinning exponent *nst* comprises contributions of red blood cell aggregation and deformability:

$$n\_{\mathfrak{sl}} = n\_{\mathfrak{sl}\mathfrak{J}} + n\_{\mathfrak{d}cf} \tag{17}$$

where each component is described by a power law:

$$m\_{\text{ag\%}/\text{def}} = \beta\_{\text{ag\%}/\text{def}} \,\gamma\_{\text{ag\%}/\text{def}}^{\prime} e^{-\nu\_{\text{ag\%}/\text{def}}}\prime,\tag{18}$$

with the empirical coefficient *β* and *ν*, and the non-dimensional shear rate *γ* = 1 + (*λγ*˙)*νg* , as defined by Carreau and Yasuda [31], with a time constant *λ*, and *νg* = 2. This formulation ensures finite *nst* at zero shear.

In the 5-component form the aggregation and deformation influences on the shearthinning exponent are combined into a single power law, due to the limited data on these effects:

$$m\_{st} = \beta \gamma^{\prime - \nu}.\tag{19}$$

The model proposed by Hund et al. allows for inclusion of the influence of large molecule concentration (proteins polysacharides, lipids), as well as fibrinogen, and temperature on the constitutive model. Due to a lack of data these are not included in the 5-parameter model.

#### 2.2.4. Yeleswarapu-Wu Model

This model is based on a visco-elastic Oldroyd-B model developed by Yeleswarapu et al. [27,32]. In this study, the visco-elastic effects are neglected, only the shear-thinning behaviour and haematocrit dependency are implemented. The shearthinning behaviour follows a modified Carraeu-type model based on a mixture model by Jung et al. [33].

The model is based on a mixture model and thus the viscosity is decribed as a function of plasma viscosity *μP* and red blood cell viscosity *μrbc* [32]:

$$
\mu\_{\rm mix} = (1 - \phi)\mu\_P + \phi\mu\_{\rm rbc\,\prime} \tag{20}
$$

where the red blood cell viscosity is described as:

$$
\mu\_{\infty}(\phi) + \left(\mu\_0(\phi) - \mu\_{\infty}(\phi)\right) \frac{1 + \ln(1 + k\dot{\gamma})}{1 + k\dot{\gamma}},\tag{21}
$$

where, in this implementation, *k* is a constant model parameter, and *μ*0 and *μ*∞ are modelled as third order polynomials of *φ*:

$$
\mu\_0 \quad = \quad a\_1 \phi + a\_2 \phi^2 + a\_3 \phi^3 \tag{22}
$$

$$\mu\_{\infty} \quad = \quad b\_1 \phi + b\_2 \phi^2 + b\_3 \phi^3 \tag{23}$$

#### 2.2.5. Casson-Merrill Model

The Casson model [34] is a classical non-Newtonian model in which the viscosity is modelled as:

$$
\mu = \left(\sqrt{\mu\_{\infty}} + \sqrt{\frac{\tau\_0}{\dot{\gamma}}}\right)^2,\tag{24}
$$

where *μ*∞ is the Casson viscosity (asymptote at high shear rate) and *τ*0 is the yield stress. The yield effect means that this model has a singularity at zero shear, leading to infinite viscosity. While there is an argument that blood does exhibit yield at slow time scales and low shear, this effect will typically make this type of model unsuited for numerical simulation within a generalised Newtonian approach with a local effective viscosity due to numerical instability.

For blood, Merill et al. gave the expressions for *μ*∞ and *τ*0 as [24,26]

$$
\mu\_{\infty} = \quad \left(\frac{\mu\_{pl}}{(1-\phi)^a}\right) \tag{25}
$$

$$\pi\_0 = \left. \beta^2 \left[ \left( \frac{1}{1-\phi} \right)^{a/2} - 1 \right]^2 \right. \tag{26}$$

with the fitting parameters *α* and *β*.

#### 2.2.6. Characteristics of Rheology Models

All viscosity model parameters were fitted to experimental data for varying levels of haematocrit in ADC plasma reported by Brooks [28]. While these data are for steady state shear only, it is still considered on of the best datasets for blood rheology data and is used in the majority of work on blood rheology. The parameters were fitted using a Levenberg–Marquardt least squares fit, implemented in Scientific Python (SciPy), using the MINPACK library. Table 1 shows the parameter sets for the different models, Figure 1 shows the comparison of model results and experimental data. All models show a good fit to the experimental data in the range were experimental data is available (*γ*˙ > 0.15 s−1), while the behaviour for low shear stress varies between the models. The Casson model shows a singularity for zero shear (yield stress behaviour), while the other models all have finite viscosity for zero shear. However, the values at low shear vary widely. For *γ*˙ = 10−2 s−1, the range of relative viscosity is between *μ*/*μP* = 71.4 to 936. This variation will heavily influence the behaviour at low shear rate, e.g., on the axis of the flow.

Figure 2 shows the parameter fit for the Quemada model, where the classical Cokelet fit [30] exhibits singularities at 12.2%, 18%, 73.1%, and 85.6% for zero shear. The Das variation [24] improves on this, but the original parameter set by Das still shows a singularity for 80.4% haematocrit. The new parameter fit performed in this study removes

the singularities completely and shows monotonous behaviour for the whole range of haematocrit concentrations and shear rates.

**Figure 1.** Comparison of non-Newtonian rheology models: Quemada (**a**), modified 5-parameter Krieger (**b**), Yeleswarapu-Wu (**c**), and Casson-Merril (**d**) model. All model parameters have been fitted to Brooks' data. Dots: experimental data (Brooks), dotted lines: model equations, data series: haematocrit concentration.

**Table 1.** Viscosity model parameters. Levenberg–Marquardt least squares fit (SciPy, MINPACK), to Brooks's data (all viscosities calculated in Pa s), *μP* = 1.23 × 10−3 Pa s.

**Figure 2.** Comparison of Quemada parameterisation for zero (**a**) and low (0.15 s−1) (**b**) shear rate. The classic Cokelet parameter set shows singularities at 12.2%, 18%, 73.1%, 85.6%, the modified parameterisation by Das improves on this, but still shows a singularity for 80.4% haematocrit. The current parameter set removes the singularity and shows monotonous behaviour.

#### *2.3. Implementation*

The fundamental equations for mass and momentum conservation were implemented using the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) [35] method for steady state, and the PISO (Pressure-Implicit with Splitting of Operators) [36] and PIMPLE (combining PISO and SIMPLE) methods for transient simulations.

Discretision is typically second order in space and time. The code supports all discretisation methods that are supported in the FOAM library.

The haematocrit transport Equation (5) is implemented as a scalar transport Equation (Listing 1), solved outside of the SIMPLE loop. The Laplacians in *φ* are implemented implicitely (*fvm*) as diffusion terms, while the source terms in *γ*˙ and *μ* are calculated explicitely (*fvc*).

**Listing 1.** Haematocrit transport equation as implemented in Openfoam.

```
gammaDot = pow(2,0.5)*mag(symm(fvc::grad(U)));
sourceC = fvc::laplacian(Kc*sqr(a)*sqr(H), gammaDot);
sourceV = fvc::laplacian(Kmu*sqr(a)
                           * gammaDot*sqr(H)/laminarTransport.nu(),
                           laminarTransport.nu());
fvScalarMatrix HEqn
                (
                fvm::ddt(H)
                + fvm::div(phi, H)
                - fvm::laplacian(Kc*sqr(a)*H*gammaDot, H)
                ==
                sourceC
                + sourceCnonlin
                + sourceV
                );
```
For steady state (SIMPLE) and transient cases with the PIMPLE algorithm, underrelaxation is required. Typically the underrelaxation factor that is required can be estimated from the order of magnitude of the ratio between collision radius and vessel radius. Stable simulation has been achieved for relaxation factors of 0.1 log(*O*(*a*/*R*)), e.g., a radius *R* = 50 μm and collision radius of 3.5 μm will require an underrelaxation factor of ≈0.1 with no underrelaxation for the final iteration. The PISO algorithm does not use underrelaxation and requires a time estimated as the smaller of (a) timestep estimated from the Courant number (*Co* < 1), and (b) time step calculated based on a Courant number scaled with the migration velocity (instead of the convective velocity).

The discretisation schemes used in the calculations presented in this paper are: second order Euler backward in time and second order (Gauss linear, and Gauss linear upwind for advective terms) in space, gradients are approximated using the least squares theme.

Rheology models are implemented as quasi-Newtonian, with calculation of local cell viscosity based on the shear rate and haematocrit value in the cell from the previous iteration/time step. The new rheology models that are implemented at the time of writing are the standard Krieger–Dougherty, the modified 5-parameter Krieger, the Yeleswarapu– Wu, and the Quemada model.

#### **3. Results**

All results shown in this paper are for fully developed pipe flow, with periodic boundary conditions between outlet and inlet, with prescribed average velocity. The radius of the pipe varies between 50 μm and 5 mm, to represent typical vessel diameters. The pipe length is two diameters.

#### *3.1. Verification and Influence of Mesh Type*

The verification case for the implementation is a pipe of radius 50 μm, average velocity *V* = 0.0065 m s−1. The rheology model used in the verification case is the Krieger– Dougherty model to allow comparison to the analytical solution [14] (no analytical solution available for the non-linear terms in the shear-stress and concentration dependent models). Model parameters for the Krieger–Dougherty model are *Kc* = 0.41, *Kμ* = 0.62, *φ*∗ = 0.68, and *n* = 1.82.

The simulation was performed for different meshes, Figure 3, (a) an axisymmetric (2D) wedge with 50 cells resolution in radial direction, (b) a hexahedral, block structured mesh— 50 cells radial, and (c) a polyhedral mesh with boundary layer inflation with ~60 cells across the diameter—this type of mesh is common in the simulation of vascular flow in patient specific geometries. The given resolutions were chosen based on a mesh convergence study and realistic mesh resolutions as typically used in vascular simulations. The migration model requires a mesh that is of similar resolution as meshes that aim at resolving wall shear stress (WSS) and WSS derived metrics.

**Figure 3.** Mesh topology for the verification of the model: axisymmetric wedge, 50 cells radial; hexahedral, block-structured, 50 cells radial; polyhedral with boundary layer extrusion, 60 cells diameter.

Figure 4 shows the results for the different meshes in comparison to the analytical solution of the migration model with the Krieger–Dougherty model. The axisymmetric two-dimensional and the hexahedral three-dimensional meshes show excellent agreement, with only a slight rounding of the peaked analytical solution at the axis. The polyhedral three-dimensional mesh also shows good agreement, but the additional numerical diffusion blunts the profile at the axis, the concentration close to the wall is well represented.

**Figure 4.** Steady state particle distribution (**a**) and velocity profiles (**b**) for different mesh types, compared with analytical solution for particle distribution by Krieger et al. Parameters: fully developed pipe flow, *R* = 50 μm, *V* = 0.0065 m s−1, *Kc* = 0.41, *Kμ* = 0.62, *n* = 1.82, *φ*∗ = 0.68, Standard Krieger–Dougherty Model.

#### *3.2. Length and Time Scale Dependency*

#### 3.2.1. Wall Shear Strain Scaling

The parabolic velocity profile for a Newtonian flow is given as:

$$v = -2\, V \left(\frac{r^2}{R^2} - 1\right),\tag{27}$$

where *V* is the average velocity.

Therefore, the velocity gradient in radial direction is:

$$\frac{\partial v}{\partial r} = -\frac{4\,\mathrm{V}r}{R^2}.\tag{28}$$

So the gradient at the wall (*r* = *R*) scales with *V* and *R*−1. The velocity is, therefore, scaled with *R*, such that the wall velocity gradient is constant. The Reynolds number scales with *R*2. For the given values of *R* = 0.05, 0.5, 5 mm, *V* = 0.0065, 0.065, 0.65 m s−1, the wall velocity gradient is constant at *γ*˙ *w* ≈ 650 s−1, to cover the significant three decades of shear strain magnitude for shear-thinning non-Newtonian blood models.

The steady state particle distribution profile is independent of the length scale and the diameter ratio. It will only depend on the ratio of *Kc*/*Kμ*. Figure 5 shows steady state profiles for a range of diameters from 0.1 to 10 mm. The computational effort for the particle migration model, however, scales with *R*2/*a*2, with *R*, the vessel radius, and *a*, the particle collision radius. While the small diameter *D* = 0.1 mm case is fully converged after around 104 iterations, the *D* = 10 mm case requires 106 iterations. This corresponds to the diffusion timescales.

**Figure 5.** Steady state particle distribution (**a**) and velocity profiles (**b**) for different diameters. Parameters: fully developed pipe flow, *R* = 0.05, 0.5, 5 mm, *V* = 0.0065, 0.065, 0.65 m s−1, *Kc* = 0.41, *Kμ* = 0.62, *n* = 1.82, *φ*∗ = 0.68, Standard Krieger–Dougherty Model.

#### 3.2.2. Kinematic and Particle Migration Timescales

Blood flow with particle migration is governed by several different time scales for flow kinematics and particle migration. The timescale for the development of the velocity profile (kinematic timescale) is

$$
\pi\_k = \frac{R^2}{\nu}.\tag{29}
$$

The timescales for the development of the particle migration profile can be derived from the particle migration flux diffusion terms as:

$$
\pi\_{c\phi} = \frac{R^2}{K\_c a^2 \phi \dot{\gamma}'} \tag{30}
$$

$$
\tau\_{c\uparrow} = \frac{R^2}{K\_c a^2 \phi^2} \tag{31}
$$

$$\pi\_{\mu} = \frac{R^2}{\mathcal{K}\_{\mu}a^2 \dot{\gamma} \frac{\partial(\ln \mu)}{\partial \dot{\Phi}}}. \tag{32}$$

The kinematic timescale scales with *R*2/*ν*, while the particle migration timescales scale with the square diameter ratio *R*2/*a*2, where *R* is the pipe radius, and *a* is the particle (collision) radius.

The kinematic viscosity, *ν* ≈ 3 × 10−6 m2 s−1, while for an average collision radius of red blood cells of *a* = 3.5 μm, the particle migration diffusion coefficients are of the order of 10−9–10−11 m2 s−1. This means the particle migration happens on timescales that are three orders of magnitude greater than the kinematic timescales.

Figure 6 shows the temporal development of the particle distribution and non-Newtonian velocity profile. The flows were initialised with a fully developed parabolic velocity profile and a uniform particle distribution of *φ* = 0.45 volume fraction. The 0.1 mm case has reached steady state conditions within 0.5 s, the 1.0 mm case shows significant particle migration after physiologically relevant times, while the 10 mm case does show only minimal migration after 10 s. It can be seen that temporal scaling follows the predicted *R*2/*a*2 scaling factor.

#### *3.3. Variation of Rheology Model and Collision Parameter Ratio*

As is obvious from Equation (5), the particle migration is strongly dependent on the viscosity model and the balance between collision and viscosity driven migration, as expressed in the model parameters *Kc* and *Kμ*.

While the magnitude of *Kc* and *Kμ* controls the magnitude of the migration pressures and thus the temporal response of the system, the concentration profile only depends on the balance between collision and viscosity driven fluxes. This balance is expressed by the ratio between the parameters *Kc*/*Kμ*. Figure 7 shows the haematocrit profiles as they develop for different viscosity models—Krieger–Dougherty (K-D), Quemada (Q), Yeleswarapu–Wu (Y), modified 5-parameter Krieger model (K5), and varying K-ratios *Kc*/*Kμ* = 0.4 to 0.75.

Compared to the verification K-D case with K-ratio of 0.66, it can be seen that a shift in the balance to higher influence of the collision frequency (higher K-ratio) steepens the profile, while a lower K-ratio, i.e., a shift of the balance to the resistive influence of the viscosity increase in the low shear region causes a flatter profile.

Comparing the different viscosity models clearly shows the main difference in the core region, where the strong variation in the low shear behaviour, discussed earlier, leads to a strong variation in the relative viscosity gradient (last term in Equation (5)). It is obvious that there is a need for further study and comparison with experimental or mesoscale modelling data to find realistic parameters for each of the potential viscosity models. Especially the modified Krieger model (K5) shows a, most likely unrealistic, double-bump profile at the axis.

Based on these results, the Quemada model with a K-ratio of between 0.5 and 0.6 seems to be the most promising candidate for a semi-mechanistic rheology model for blood.

**Figure 6.** Transient particle distribution (**a**,**c**,**e**) and velocity profiles (**b**,**d**,**f**) for different diameters. Parameters: fully developed pipe flow, *R* = 0.1 (**a**,**b**), 1.0 (**c**,**d**), 10 (**e**,**f**) mm, *V* = 0.0065 (**a**,**b**), 0.065 (**c**,**d**), 0.65 (**e**,**f**) m s−1, *Kc* = 0.41, *Kμ* = 0.62, *n* = 1.82, *φ*∗ = 0.68, Standard Krieger–Dougherty Model.

**Figure 7.** Steady state particle distribution (**a**,**c**,**e**,**g**) and velocity profiles (**b**,**d**,**f**,**h**) for different viscosity models and collision parameter ratios (*Kc*/*Kμ*). Parameters: fully developed pipe flow, *R* = 0.05 mm, *V* = 0.0065 m s−1, *Kc*/*Kμ* = 0.4–0.75 (given in legends). Standard Krieger–Dougherty (**a**,**b**), Quemada (**c**,**d**), Yeleswarapu (**e**,**f**), and modified 5-parameter Krieger (**g**,**h**) model.

#### *3.4. Application to Realistic Vessel Model*

The hybrid implicit–explicit implementation presented here allows for transparent change of viscosity model and does not use any artificial stabilisation, allowing for low residual solutions on realistic complex geometries.

Figure 8 shows the time-averaged haematocrit concentration near the wall and the time-averaged WSS for a realistic human common carotid artery (CCA).

The variation of haematocrit in this case is around ±1%, which is in agreement with the expected haematocrit variation for the size (*RCCA* = 3 mm).

While the haematocrit transport is not significant in this case, it shows that the implementation is also stable for physiological vessel shapes (and the reduced grid quality that is typical for clinical applications) and for highly dynamic flow. The haematocrit-dependent viscosity models can be employed, and for larger vessels, where the haematocrit transport can be neglected, the haematocrit equation can be frozen, to reduce simulation overhead.

**Figure 8.** Time averaged haematocrit (**a**) and WSS (**b**) in a realistic human carotid, 900 k cells.

#### **4. Discussion**

While previous implementations [17,19] of this class of model are limited by the fact that the viscosity term in Equation (5) is linearised in the viscosity gradient with H, our implementation avoids this by implementing the non-linear term directly, which (i) allows the use of different viscosity models and, crucially, (ii) facilitates the functional forms of transport parameter dependence on field variables we are constructing from the operation of our explicit, meso-scale models. Furthermore, the present implementation avoids the use of artificial stabilisation terms, that lead to underestimation of RBC migration [18].

The particle migration time scales with (*a*/*R*)2, where *a* is the RBC collision radius. This means that the particle migration is most relevant for small vessels of a diameter of 1 mm or lower, where the migration occurs on physiologically relevant timescales. For

larger vessels, minor effects are caused by a synergy of particle migration and secondary flows [18] (Figure 8).

#### *Limitations of the Model and Future Work*

It has to be noted that the implementation presented here uses the magnitude of the shear in the particle flux formulation. As noted by Phillips [14], this assumes an essentially one-dimensional shear state, and an isotropic response. This limits the application of the model macroscopic to flow situations where the shear tensor is aligned with the flow and the main shear in radial direction, though it also provides a convenient initial simplification to our microscopic modelling. As with isotropic turbulence modelling the isotropic migration model will over-predict migration pressure in regions with high anisotropy, e.g., stagnation points, strong acceleration, or rotational shear. The authors are currently looking to implement an explicit formulation for a localised, anisotropic shear and migration pressure tensor, similar to approaches proposed by Miller [37] or Fang et al. [15].

Clearly, the parameters of the particle migration model we have used in this work are intentional targets for assignment using data from the operation of our explicit mesoscale simulations. Of course, for immediate impact, one might calibrate them to suitable experimental data, using established phenomenology. (Though note that the improved quantitative understanding we hope to obtain from our direct, meso-scale studies might produce the further benefit of insight into the functional form of transport parameter dependence upon field variables). While such experimental data are available, albeit scarce, for rigid particles in suspension, e.g., based on nuclear magnetic resonance measurements of particle profiles, the authors are not aware of any such data for soft vesicles, in particular RBCs. Here again, we hope eventually to use our recent meso-scale models (MCLBM) which account for the cellular-scale interactions in full detail [11] to refine and parameterise a constitutive equation for concentrated suspensions of deformable vesicles and hence to inform the present continuum model parameters. The results presented here provide us with the reassurance that, the final discussions below notwithstanding, a suitably sensitive, macro-scale model exists for this undertaking.

#### **5. Software**

The continuum-scale haemorheology framework was implemented in OpenFoam, version 1912 and 2012.

haemoFoam is a modelling framework for vascular flow simulation based on FOAM, that is intended to cater for the particular requirements of haemodynamics, in particular with respect to WSS related phenomena like atherosclerosis. At the time of writing it includes:

	- **–** Krieger Dougherty (non shear-thinning);
	- **–** Modified K-D [25] (shear-thinning);
	- **–** Quemada;
	- **–** Yeleswarapu;
	- **–** Casson–Merrill;
	- **–** Carreau model (not concentration dependent, Fluent implementation).
	- **–** TAWSS, TAWSSMag;
	- **–** OSI;

Planned future features are:

haemoFoam is open-sourced under GPL3 and will be made available to interested parties upon request.

**Author Contributions:** Conceptualization, T.S. and I.H.; methodology, T.S. and I.H.; software, T.S.; validation, T.S.; formal analysis, T.S.; writing—original draft preparation, T.S.; writing—review and editing, T.S. and I.H.; visualization, T.S.; funding acquisition, I.H. and T.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by NHS Heart of England Trust grant number REC 15/WM/0164.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not appicable.

**Data Availability Statement:** The software (haemoFoam) is freely available to interested parties on github (TS-CUBED/haemoFoam). Please contact the author for testing and developer access.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


**4. Peroxiredoxin-4**

#### **4.1 Peroxiredoxin-4 and its function**

H2O2 is produced during the disulfide bond formation in protein synthesis and from a side effect of cell biology mechanism located in the ER lumen and extracellular matrix [54]. PRDX4 is localized in the endoplasmic reticulum and highly expressed in the pancreas, liver, and heart but is low in blood and brain leukocytes [5]. PRDX4 plays a role as an antioxidant by scavenging H2O2 and is a chaperone molecule to activate ER stress pathway [55]. PRDX4 is important for cell protection by reducing H2O2 to water in order to reduce oxidative stress.

In addition, PRDX4 has a role in regulating the activation of NF-ĸB transcription factor and TNF-related apoptosis-inducing ligand (TRAIL) [56]. TRAIL is a transmembrane protein involved in apoptosis extrinsic pathway by binding to receptors DR4 or DR5. A reduction in PRDX4 expression leads to TRAIL activation and apoptosis induction, whereas an increase in PRDX4 expression inhibits the TRAIL activation and leads to cell survival [57]. The expression of PRDX4 in the cytosol increases NF-ĸB activity through phosphorylation of IκB-α [58]. Meanwhile, the expression of PRDX4 in the endoplasmic reticulum reduces NF-ĸB activity. NF-ĸB plays a role in cell proliferation and survival. A study has shown that PRDX4 is bound to the endothelium and is secreted when there is a redox change in the extracellular matrix [59]. This study is supported with the presence of high PRDX4 in the serum that indicates membrane leakage due to tissue destruction and cell apoptosis [60].

#### **4.2 Mechanism of peroxiredoxin-4 as an antioxidant**

PRDX4 triggers H2O2 signaling and protects cells from oxidative stress by oxidizing thioredoxin (Trx) proteins. Trx is a group of co-factor proteins with chaperone activity to control cell homeostasis and inhibit H2O2 signaling pathway [61]. PRDX4 requires Trx peroxidase activity as an electron donor. The oxidized Trx is recycled by Trx reductase and NADPH as the primary electron donor. The main target of H2O2 molecules is two cysteine residues on PRDX4 to form catalytic peroxide mechanisms [62]. PRDX4 has two cysteine residues known as peroxidatic cysteine residue (CysP-) and resolution cysteine residue (CysR-). These two cysteine subunits have peroxidatic properties to oxidize H2O2.

When the oxidative stress level is high, CysP- acts on peroxide to form sulfenic acid (CysP-SOH). CysP- is one of the cysteine units located near the N-terminal of PRXD4. CysP-SOH then reacts with another systemic subunit known as CysR- to form a stable disulfide bond, and water molecule is produced. Trx enzyme reduces the disulfide bond and completes the PRDX4 catalytic cycle by detoxifying the peroxide and producing an active thiol. If the oxidative stress keeps increasing, –SOH is

**109**

*Silencing of Peroxiredoxin-4 in Anticancer Activity of Gamma-Tocotrienol*

oxidized to –SO2H [63]. PRDX4 in -SO2H is reduced and re-oxidized by sulfiredoxin (Srx) through the reduction of ATP to ADP. This process is reversible [1]. However, in high oxidation state, –SO2H is continuously oxidized to –SO3H, and PRDX4 will undergo hyperoxidation. This will cause inhibition of oxidation, and PRDX4 will activate other pathways to counteract the oxidation state. In conclusion, PRDX4 antioxidant activity is inhibited in high oxidation state and indirectly activates the

Studies have shown that the ER produces higher level of H2O2 than the mitochondria [65]. The main source of H2O2 in the ER lumen is the reaction of oxidase-1 in the ER (Ero1) with disulfide isomerase protein (PDI) for disulfide bond formations during protein synthesis. The disulfide bonds are important to prevent the protein from degrading after it is excreted from the ER lumen. Disulfide bonds formation involves high oxidative processes, and the unfolded bonds on damaged proteins may cause unfolded protein response (UPR). Unfolded proteins are degraded by dissolution of disulfide bonds before the protein is bound to the cytosol membrane [66]. UPR is the action of cells to balance the increasing oxidative stress and repair the function of or degrade the unfolded protein. The UPR will activate chaperone proteins of multiple pathways to balance the oxidative stress [67]. High oxidative stress in the ER will induce UPR and activate Ero1 to oxidize PDI catalytic enzymes by using oxygen molecules as the electron receiver [68]. The UPR action will form three intermediate sensors which are ATF6, IRE1α, and PERK and at the same time increase Ero1 and PRDX4 expressions as homeo-

**5. The effect of viability, apoptosis rate, and reactive oxygen species production in HepG2 cell silenced with peroxiredoxin-4 and treated**

Gamma-tocotrienol has been shown to have high antitumor activity [12] that exhibits cell cycle arrest and apoptosis in alveolar adenocarcinoma epithelial cells, A549 [70], and colorectal cancer cell HT-29 [71]. Its action is specific and causes the cancer cells to be more susceptible and sensitive to treatment [72]. A previous study by Sazli et al. [7] reported that PRDX4 was upregulated when HepG2 cells were treated with GTT. In this study, we silenced *PRDX4* gene and treated with 70 μM of GTT to determine the role of PRDX4 in anticancer activity of GTT. Finding has shown that GTT treatment reduces HepG2-shRNA-PRDX4 cell viability significantly (*p* < 0.01) compared to its control (HepG2-shRNA-PRDX4) and also compared to HepG2 treated with GTT (*p* < 0.05) (**Table 1**). This viability assay is the first line of observation to show that GTT treatment does have an effect toward the

To investigate further the cause of reduction in cell viability, an apoptosis assay was done. The findings showed that the apoptosis rate was significantly increased (*p* < 0.05) in HepG2 treated with GTT compared to HepG2 control. Besides that, apoptosis rate increased in HepG2-shRNA-PRDX4-GTT (*p* < 0.05) compared to HepG2-shRNA-PRDX4 and HepG2-GTT (*p* < 0.05). GTT is capable to induce intrinsic and extrinsic apoptotic pathways in cancer cells such are breast cancer cell line, MCF-7, and MDA-MB-231 by increasing biomarker stress in the endoplasmic reticulum [73]. It is suggested that the main reason for the apoptosis in GTTtreated group is high level of ROS production. GTT has triggered ROS production through PRDX4 activity. Studies reported that one of the anticancer activities of GTT is apoptosis induction through increases of ROS production [74]. We have postulated that GTT might trigger the ROS production through PRDX4 activity.

*DOI: http://dx.doi.org/10.5772/intechopen.88813*

ER-pressure pathway [64].

static chaperones [69].

**with GTT**

silenced *PRDX4* gene.

*Synthetic Biology - New Interdisciplinary Science*

and protein expression of targeted gene.

**4.1 Peroxiredoxin-4 and its function**

by reducing H2O2 to water in order to reduce oxidative stress.

**4.2 Mechanism of peroxiredoxin-4 as an antioxidant**

subunits have peroxidatic properties to oxidize H2O2.

**4. Peroxiredoxin-4**

cell apoptosis [60].

the plasmid carried by the lentivirus has no effect on the cell. Although shRNA functions are at the mRNA level, only changes in protein expression may contribute to phenotypic changes. Therefore, validation of the mRNA and protein expression of the targeted proteins should be carried out to ensure that the gene expression is suppressed [53]. In addition, it is essential to validate the silencing effect by using more than two sequences of targeted gene shRNA. Every segment of gene sequence has a different effect in silencing the gene. The best sequence is selected based on the higher silencing rates shown through reduction of mRNA

H2O2 is produced during the disulfide bond formation in protein synthesis and from a side effect of cell biology mechanism located in the ER lumen and extracellular matrix [54]. PRDX4 is localized in the endoplasmic reticulum and highly expressed in the pancreas, liver, and heart but is low in blood and brain leukocytes [5]. PRDX4 plays a role as an antioxidant by scavenging H2O2 and is a chaperone molecule to activate ER stress pathway [55]. PRDX4 is important for cell protection

In addition, PRDX4 has a role in regulating the activation of NF-ĸB transcription factor and TNF-related apoptosis-inducing ligand (TRAIL) [56]. TRAIL is a transmembrane protein involved in apoptosis extrinsic pathway by binding to receptors DR4 or DR5. A reduction in PRDX4 expression leads to TRAIL activation and apoptosis induction, whereas an increase in PRDX4 expression inhibits the TRAIL activation and leads to cell survival [57]. The expression of PRDX4 in the cytosol increases NF-ĸB activity through phosphorylation of IκB-α [58]. Meanwhile, the expression of PRDX4 in the endoplasmic reticulum reduces NF-ĸB activity. NF-ĸB plays a role in cell proliferation and survival. A study has shown that PRDX4 is bound to the endothelium and is secreted when there is a redox change in the extracellular matrix [59]. This study is supported with the presence of high PRDX4 in the serum that indicates membrane leakage due to tissue destruction and

PRDX4 triggers H2O2 signaling and protects cells from oxidative stress by oxidizing thioredoxin (Trx) proteins. Trx is a group of co-factor proteins with chaperone activity to control cell homeostasis and inhibit H2O2 signaling pathway [61]. PRDX4 requires Trx peroxidase activity as an electron donor. The oxidized Trx is recycled by Trx reductase and NADPH as the primary electron donor. The main target of H2O2 molecules is two cysteine residues on PRDX4 to form catalytic peroxide mechanisms [62]. PRDX4 has two cysteine residues known as peroxidatic cysteine residue (CysP-) and resolution cysteine residue (CysR-). These two cysteine

When the oxidative stress level is high, CysP- acts on peroxide to form sulfenic acid (CysP-SOH). CysP- is one of the cysteine units located near the N-terminal of PRXD4. CysP-SOH then reacts with another systemic subunit known as CysR- to form a stable disulfide bond, and water molecule is produced. Trx enzyme reduces the disulfide bond and completes the PRDX4 catalytic cycle by detoxifying the peroxide and producing an active thiol. If the oxidative stress keeps increasing, –SOH is

**108**

oxidized to –SO2H [63]. PRDX4 in -SO2H is reduced and re-oxidized by sulfiredoxin (Srx) through the reduction of ATP to ADP. This process is reversible [1]. However, in high oxidation state, –SO2H is continuously oxidized to –SO3H, and PRDX4 will undergo hyperoxidation. This will cause inhibition of oxidation, and PRDX4 will activate other pathways to counteract the oxidation state. In conclusion, PRDX4 antioxidant activity is inhibited in high oxidation state and indirectly activates the ER-pressure pathway [64].

Studies have shown that the ER produces higher level of H2O2 than the mitochondria [65]. The main source of H2O2 in the ER lumen is the reaction of oxidase-1 in the ER (Ero1) with disulfide isomerase protein (PDI) for disulfide bond formations during protein synthesis. The disulfide bonds are important to prevent the protein from degrading after it is excreted from the ER lumen. Disulfide bonds formation involves high oxidative processes, and the unfolded bonds on damaged proteins may cause unfolded protein response (UPR). Unfolded proteins are degraded by dissolution of disulfide bonds before the protein is bound to the cytosol membrane [66]. UPR is the action of cells to balance the increasing oxidative stress and repair the function of or degrade the unfolded protein. The UPR will activate chaperone proteins of multiple pathways to balance the oxidative stress [67]. High oxidative stress in the ER will induce UPR and activate Ero1 to oxidize PDI catalytic enzymes by using oxygen molecules as the electron receiver [68]. The UPR action will form three intermediate sensors which are ATF6, IRE1α, and PERK and at the same time increase Ero1 and PRDX4 expressions as homeostatic chaperones [69].


**1. Introduction**

Trees, including those outside of forests (i.e., orchards, gardens, trees in urban environments), are terrestrial plants that provide several ecosystem services essential to life and the economy. Trees contribute to carbon sequestration, animal and insect biodiversity maintenance, lignocellulosic biomass production for the industry, food supply (i.e., fruits and seeds), tourism (i.e., forest nature trails), air filtration in urban areas, and wood production for construction. Given their importance and the benefits trees provide, understanding the factors related to their health is pivotal for proper human intervention in forests, woodlands, and orchards.

Tree health depends on numerous factors, and it is now understood that their health, as well as other multicellular organisms, is closely related to the structure of the microbial community of the holobiont (the host and its associated microorganisms) [1]. The theory of plant holobiont seeks to shed more light on the plant as a meta-organism,

in which plant growth, health, and productivity are closely related to the composition and functions of the microbiota inhabiting different niches in the plant [2]. Within this perspective, trees are not autonomous entities, and their health can be considered a function of the plant microbiota (set of all associated microorganisms) or the plant microbiome (set of all associated microbial genomes) [3–4].

Trees establish different symbiotic relationships with a plethora of microorganisms. Therefore, the plant microbiome, also known as a phytobiome, is composed of commensal, mutualistic, and pathogenic microorganisms [5]. The diverse microbial composition contributes to host-microbial homeostasis, meaning the microbial community structure can modulate and maintain tree health over time. It is now known that the colonization of plant tissues by a pathogenic microorganism and disease development may be related to an imbalance of host-microbial homeostasis. In some cases, this imbalance results from atypical environmental conditions or the loss of partner microbes or ancestor microorganisms that are key to tree health.

The microorganisms that are part of the phytobiome include a wide diversity of prokaryotes distributed in different phyla (Proteobacteria, Actinobacteria, Firmicutes, Bacteroidetes, etc.), in addition to fungi (including mycorrhizae), protozoa, nematodes, and viruses. In general, the rhizospheric soil and endosphere are the environments most densely populated by microorganisms; nevertheless, the phylloplane is also an important environment for tree-associated microbial communities. Partner microorganisms can contribute directly and indirectly to tree health through nutrient provision (i.e., phosphorus solubilization and nitrogen fixation), plant hormone production (i.e., gibberellins, auxins), pathogen inhibition via competition, and the inhibition of systemic host resistance.


**4.4. Risk constraint formulation**

As already mentioned before, the financial decision process can be modeled by means of a MDP. Naming profit *B w* for each sample the sum of income, cost and transaction cost over all instruments, the objective function (56) becomes:

$$\max\_{\mathbf{x}} \left[ \mathbb{E} \left[ \sum\_{t} \left[ \boldsymbol{\mathcal{B}}^{\mathrm{w}} \left( \mathbf{x}\_{t}, \mathbf{x}\_{t-1} \right) \right] \right] \right] \tag{70}$$

Exchanging the order of the expectation operator and the summation, and expanding the summation:

$$\max\_{\mathbf{x}} \left[ \mathbb{E} \left[ \mathcal{B}\_1^{\mathbf{w}} \left( \mathbf{x}\_1, \mathbf{x}\_0 \right) \right] + \sum\_{t=2}^n \mathbb{E} \left[ \mathcal{B}\_t^{\mathbf{w}} \left( \mathbf{x}\_t, \mathbf{x}\_{t-1} \right) \right] \right] \tag{71}$$

where *n* is the total number of time periods considered.

At this point and considering the first term depends only on the initial market position *x*0 and the first strategy decision *x*1, the maximization can be decomposed using the Bellman´s Principle of optimality as follows:

Risk-Constrained Forward Trading Optimization by Stochastic Approximate Dynamic Programming http://dx.doi.org/10.5772/57466 115

$$\max\_{\mathbf{x}\_{l}} \left[ \mathsf{E} \left[ \boldsymbol{\mathcal{B}}\_{l}^{\mathrm{w}} \left( \mathbf{x}\_{l}, \mathbf{x}\_{0} \right) \right] + \max\_{\mathbf{x}\_{t:2...n}} \sum\_{t=2}^{n} \mathsf{E} \left[ \boldsymbol{\mathcal{B}}\_{t}^{\mathrm{w}} \left( \mathbf{x}\_{t}, \mathbf{x}\_{t-1} \right) \right] \right] \tag{72}$$

If *x* is independent of the Monte Carlo sample, the terms inside the summation over all future periods (*t* =2…*n*) are simply the expected profit for the period *t* after a decision *xt*−1→ *xt*. However, the model should take into account that future strategy decisions may be different for each Monte Carlo sample, accounting for adjustments the decision-maker almost certainly would execute to face specific scenarios. Then, the expected profit *Bt* ¯ and the decision itself for future decision stages will depend on a set of variables *st*, which represent the variables considered by the decision maker in order to adapt the strategy to a particular situation and the equation (72) becomes:

$$\max\_{\mathbf{x}\_{l}} \left[ \overline{B\_{l}} \left( \mathbf{x}\_{l}, \mathbf{x}\_{0}, s\_{l} \right) + \max\_{\mathbf{x}\_{r=2...n}} \sum\_{t=2}^{n} \overline{B\_{t}} \left( \mathbf{x}\_{t}, \mathbf{x}\_{t-1}, s\_{t} \right) \right] \tag{73}$$

Defining the continuation value functions *Vt* as:

is below the prevailing spot prices, following the decision model of equation (42). Variable costs of generation are assumed linear with power output, i.e. marginal costs are constant.

The operation-failure cycles of the generating unit are obtained from a chronological Marko‐ vian stochastic simulation. For each spot price sample, a time series of power output is synthesized for every generation unit. The hourly power output is simulated o following three

**1.** Based on failure and repair rates defined by the state the unit resided in the previous hour, a random failure is simulated [12]. If a failure is in place, the output power is set to zero

**2.** The dispatch of the unit is simulated, taking into account the marginal cost of generation and the prevailing sample spot price at that time interval. Here perfect foresight of the spot price is assumed in order to decide the dispatch and fulfill the minimal generation

**3.** If dispatched, other unit's technical restrictions are fulfilled, e.g. ramping capabilities.

depending on whether the unit is generating or is in stand-by.

114 Dynamic Programming and Bayesian Inference, Concepts and Applications

instruments, the objective function (56) becomes:

where *n* is the total number of time periods considered.

Principle of optimality as follows:

This chronological stochastic model reproduces with accuracy the dynamics involved in failure and repair cycles of generators, giving the possibility to select different failure rates

As already mentioned before, the financial decision process can be modeled by means of a MDP. Naming profit *B w* for each sample the sum of income, cost and transaction cost over all

Exchanging the order of the expectation operator and the summation, and expanding the

1 10 ( ) ( ) 1 2

At this point and considering the first term depends only on the initial market position *x*0 and the first strategy decision *x*1, the maximization can be decomposed using the Bellman´s

=

é ù ê ú é ùé ù + ë ûë û ë û å *n*

E *B xx* (70)

E E *B x x B x x* (71)

max ( ) 1 , é ù é ù ê ú ê ú é ù ë û ë û ë û å *w t t x t*

max , , -

*w wt tt x t*

steps:

times.

summation:

for this particular hour.

**4.4. Risk constraint formulation**

$$V\_t\left(\mathbf{x}\_t, \mathbf{x}\_{t-1}, \mathbf{s}\_t\right) = \overline{B\_t}\left(\mathbf{x}\_t, \mathbf{x}\_{t-1}, \mathbf{s}\_t\right) + \max\_{\mathbf{x}\_{t+1}} \left[V\_{t+1}\left(\mathbf{x}\_t, \mathbf{x}\_{t+1}, \mathbf{s}\_{t+1}\right)\right] \tag{74}$$

The maximization can be solved by a set of recursive maximizations, each one solving only one decision stage:

$$\max\_{\mathbf{x}\_t} \left[ \ \ \ \mathbf{V}\_t \left( \mathbf{x}\_t, \mathbf{x}\_{t-1}, \mathbf{s}\_t \right) \right] \tag{75}$$

With this model, the optimization can be decomposed in steps and the dynamic nature of a strategy can be accurately replicated. Despite the fact future trading decisions *xt*=2…*n* are considered and optimized, the practical product of this procedure is the new optimal reba‐ lanced state *xt*=1 starting from the previous trading position *xt*=0. The further trading positions are only optimal given the current information available and should be reconsidered later. Therefore, each new trading position (*xt*=2, *xt*=3, …, *xt*=*n*) should be the product of a similar optimization incorporating the additional market information available immediately before.

The value functions provide the expected continuation value within a state space defined by the state variables, *xt*, *xt*−1 and *st* 3 . However, the continuation functions, which are essential to solve the optimization problem, are unknown beforehand. It is here that the ADP approach is introduced to approximate the value and risk functions for the state space.

3 There are several other decomposition methods, some of which exclude the decision as state space variable defining the value functions in a post-decision state space. These approaches make the step maximization sometimes harder but present advantages such as a state space of fewer dimensions. See for example [2].

The set of constraints remains the same as they were already defined for each period *t*. Nevertheless, special considerations should be made to calculate de risk and to fulfill the riskconstrained optimization. As *Vt* only account for the expected profit, risk functions *Rt* should also similarly be calculated for the same state space in order to enforce the risk constraint.

be enforced at the early iterations of the ADP algorithm, simply because the sample size of the

Risk-Constrained Forward Trading Optimization by Stochastic Approximate Dynamic Programming

http://dx.doi.org/10.5772/57466

117

With the objective of validating the results of the proposed ADP algorithm, a first simple exemplary case is considered in which a thermal generator sells energy in the spot market and in a future contract. The results of the ADP algorithm were compared with the results of a

The fractions of energy sold in the spot market and in a quarter future contract were optimized considering that the future can be traded during the delivery period. The optimization determines three decision stages during this period, one at the beginning of each month, consisting on sell or buy energy in the future market based on the previous state. The spacestate previous to each decision is defined only by the level of future already sold, in order to

The results obtained by solving the problem by means of the ADP and DP algorithms are presented in Figure 7. The plots represent the expected profit and the downside risk measured by the CVaR of the optimal strategy as a function of the initial state, i.e. the energy already committed in the future contract at the initial stage. An excellent agreement between the optimal strategies obtained by ADP and DP is evidenced, validating the proposed approach.

It can be noticed that the expected profit rises as the amount of energy sold forward increases. This is caused by the risk premium paid to the generator in the future market, i.e. the mean future prices are higher than the mean spot prices. Additionally, the transaction costs are not compensated by the risk premium; therefore the best trading strategy is to maintain similar involvement in the future market to the initial level without rebalancing the portfolio. This is illustrated in Figure 8, where the optimal decisions for the first month are practically the same

It is noteworthy to observe in Figure 7 that financial risk lessens when forward contracting in the future market increases. This means that for the conventional generator considered, which present a high availability, the delivery risk is lower than the risk of not being dispatched in the spot market. The behavior of the risk curve is closely related to the unit's failure and reparation rates and to the marginal production costs. Generators with low marginal costs are in the first places of the dispatch merit order, and hence the risk of not being dispatched is low. Moreover, high failure rates imply also a higher delivery risk. Out of these relations arise a broad number of risk curves that differ from one generator to another and suggest that considerable risk mitigation by aggregating different generators in a portfolio is possible.

when solving with a conventional DP and an ADP approach.

conventional DP algorithm for which the space-state was discretized appropriately.

dataset is small to get a consistent and statistically converged risk value.

**5. Numerical case study**

keep tractable the DP problem.

**5.2. Validation results**

**5.1. Algorithm validation**

By definition, the linear regression will deliver an approximation that minimizes the mean square error on the entire dataset. In a stochastic setting where the same inputs leads to several different outputs, the regression will accurately estimate the expected value for a given set of inputs, provided the sample size and the approximation order are appropriate. This fits perfect for the case of the value function but for the approximation of the risk function some problems arise.

Let suppose that the CVaR is chosen as risk metric. As the algorithm progresses, new data points are collected, i.e. a set of state variables and its corresponding simulated profits for the period. For approximating the CVaR associated to a particular point in the state space, one approach is to select a subset from the dataset whose input variables are "close" to the point. Then the CVaR is calculated first by sorting the profit values and then taking the mean of those below the specified α-quantile. Now, let suppose that a new data point is simulated and an update of the CVaR approximation is needed. To do so, the process described must be repeated, but now including the new data point. This simple approach has large disadvantages: all the data points must be stored and the mean is not easily updated as old data may be excluded or included of the zone below the α-quantile. These drawbacks are caused by the fact that the CVaR is quantile-based. To solve these difficulties, another solution is envisioned. Instead of approximating directly the CVaR, another risk measure is used to approximate the CVaR within the space state. The risk measure used is called Relative Lower Semideviation (RLS) and it is moment based instead. Hence, it can be updated more easily and without needing to store the entire dataset. These types of risk measures are described in detail in [13-16] and for a stochastic profit *Pt* have the form:

$$RLS\left(B\_t\right) = -\mathbb{E}\left[\left.B\_t\right] + a \cdot \sigma\_p\left^-\left[\left.B\_t\right]\right\|\right.\tag{76}$$

$$\sigma\_p \, ^- \begin{bmatrix} B\_t \end{bmatrix} = \left( \mathbb{E} \left[ \max \left( \mathbb{E} \left[ B\_t \right] - B\_t, 0 \right)^p \right] \right) \bigvee\_p \tag{77}$$

where the equation (77) is the negative semideviation of degree *p* of the stochastic profit *Pt*.

It can be proven that these moment-based risk measures are coherent if 0≤*a* ≤1 and *p* ≥1. To approximate a CVaR with a 5%-quantile, in this work the parameters used are *a* =1 and *p* =9.5. To compute the approximation, two linear regressions were used to calculate the expectations on the profit and on the negative deviations, which can be updated using the same method proposed for the value functions. In a Monte Carlo scheme, a large amount of data is needed to compute reliable risk estimations. Therefore, the risk constraint should not be enforced at the early iterations of the ADP algorithm, simply because the sample size of the dataset is small to get a consistent and statistically converged risk value.


**Author details**

Etefia Etefia Department of Medical Laboratory Science, University of Calabar, Calabar, Nigeria

\*Address all correspondence to<EMAIL_ADDRESS>
© 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.