Open access peer-reviewed chapter

# Nonlinear Oxygen Transport with Poiseuille Hemodynamic Flow in a Micro-Channel

Written By

Terry E. Moschandreou and Keith C. Afas

Reviewed: November 20th, 2019 Published: December 26th, 2019

DOI: 10.5772/intechopen.90575

From the Edited Volume

## Nonlinear Systems

Edited by Walter Legnani and Terry E. Moschandreou

Chapter metrics overview

View Full Metrics

## Abstract

In a recent paper by the authors, a well-known governing nonlinear PDE used to model oxygen transport was formulated in a generalized coordinate system where the Laplacian was expressed in metric tensor form. A reduction of the PDE to a simpler problem, subject to specific integrability conditions, was shown, and in the present work, a novel approximate analytical solution is obtained in terms of the degenerate Weierstrass P function using a compatibility relation through the factorization of the reduced almost linear ode and subject to similar boundary conditions for a microfluidic channel used in recent work by the authors. A specific form of the initial equation which was reduced has been used by Nair and coworkers describing the intraluminal problem of oxygen transport in large capillaries or arterioles and more recent work by the corresponding author describing the release of adenosine triphosphate (ATP) in micro-channels. In the present problem, a channel with a central core, rich in red blood cells, and with a thin plasma region near the boundary wall, free of RBCs is considered.

### Keywords

• almost linear ODE
• Poiseuille flow
• oxygen transport
• Painlevé analysis

## 1. Introduction

Various biophysical phenomena are modeled using nonlinear differential equations. Such is the case of a model used by Nair et al. [1, 2, 3] to describe oxygen transport in large capillaries [1, 2, 3]. This model incorporates two regions of blood flow. One is a core region of the blood with RBCs present, and in this core, oxygen dissociates into blood oxyhemoglobin. The velocity of blood in the core region is a function of the plasma velocity and rate of oxygen dissociating. The second region is a thin strip of flowing plasma with no RBCs at the wall of the micro-fluidic channel. The appropriate Robin condition and no-flux conditions are incorporated at and around a permeable membrane with oxygen transport through the membrane. Since there are two distinct regions of flow of liquid, RBCs with plasma and plasma alone, it is necessary to match the rate of change of partial pressure of oxygen at the common boundary of the liquid in each of the two regions. This kind of model has been used previously by Moschandreou et al. [4] studying the influence of tissue metabolism and capillary oxygen supply on arteriolar oxygen transport. In that study a numerical approach was used to solve the governing equations. In the present work we seek an analytical solution for a different highly nonlinear problem in a micro-fluidic channel. Flows that are fully developed at inlet were studied by Ng [5] who studied oscillatory dispersion in a tube with chemical species undergoing linear reversible and irreversible reactions at the tube wall. Relative importance to the present work is that fully developed flow occurs at inlet of channel with boundary conditions specified on the wall. No inlet mass transport is specified at inlet of channel similar to [5] but unlike [6], for example. Using Painlevé analysis [7], certain nonlinear second-order ordinary differential equations (ODE) can be factorized and solved. We consider the model of Nair et al. [1, 2, 3], where the nonlinear PDE reduced in [8], to an “almost linear” second-order ode is considered. It is well-known that the Weierstrass P function in its series form is problematic to use in computational work due to very slow convergence of numerical methods. It is the aim of the present work to show that a degenerate form of the special function can be used as in [9] in the reduction of the nonlinear PDE in [1, 2, 3]. A thorough and recent review of oxygen control with microfluidics has been carried out in [10] and all of its references within. In this work we see how the microscale can be leveraged for oxygen control of RBCs.

## 2. General tensorial mass transport

Regardless of the kinematics of a surface (dynamic or stationary), all surfaces, S=∂Ω, enclosing a solid volume, Ω, obey the following intuitive conservation relation for an enclosed observable mass, mo of some arbitrary substance:

ddtmo+SjodS=Σ,E1

where jo is the flux of the observable out or into the surface and Σ represents the net increase or decrease in the observable’s mass.

The relation states intuitively that any change of the observable’s mass within the solid, plus all observable mass entering or leaving the boundary, should represent the net change in the mass of the object.

Any mass transport can be derived from the above relation converted into the differential form. We first recognize that the observable’s mass can be represented through the observable’s density:

mo=ΩρodΩ,

In addition, we can make the same statement about the net equilibrium constant. Suppose there is a local equilibrium density, σ, such that

Σ=ΩσdΩ.

We then can obtain a full integral form of the conservation relation:

ddtΩρodΩ+SjodSΩσdΩ=0.E2

In general, it can be shown that for a general surface that is moving, the time derivative of a volume integral defined over the dynamic volume enclosed by the surface can be summarized as [11]

ddtΩψdΩ=ΩtψdΩ+SC˜ψdS,E3

where C˜ is defined as the normal projection of the surface’s perpendicular speed, also named, the surface velocity. This is encapsulated (using Einstein summation convention) by the tensorial equation:

C˜=VN=ViNi,

where N=ZiNi is the unit normal to the surface and Zi is the contravariant basis for the coordinate space. Thus, we can simplify our conservation relation:

ΩtρodΩ+SC˜ρodS+SjodSΩσdΩ=0,

We can also simplify the vector surface element, dS=NdS, and convert the flux term into a tensorial formation, by recognizing jo=joiZi:

ΩtρodΩ+SC˜ρodS+SjoiNidSΩσdΩ=0.

Finally, we will use the definition of the surface velocity and combine the two surface integral terms into one:

ΩtρodΩ+SViρo+joiNidSΩσdΩ=0.

We can use Gauss’ divergence theorem on the surface integral term and unite all the terms under one volume integral:

Ωtρo+iViρo+joiσdΩ=0.

Using the localization theorem, the integrand must be zero inside the volume integral, and we obtain the differential form of the conservation relation:

tρo+iViρo+joiσ=0.E4

We can simplify the equation, further by considering a particular form of the flux. In this, we consider advective flux (flux due to bulk movement of an observable’s mass) and diffusive flux (flux due to a concentration gradient). Using advective formulas and Fick’s first law, we obtain the flux to be

joi=voiρoDoiρo,E5

where vo=voiZi is the velocity of the observable within the volume. Substituting these into the differential conservation relation, we obtain

tρo+iViρo+voiρoDoiρoσ=0.

We simplify the equation, expanding the covariant derivative to obtain the final form of the conservation relation:

tρo+iVi+voiρo=iDoiρo+σ.E6

This form can also be put into an invariant tensorial form by utilizing the invariant time derivative operator ̇ from the calculus of moving surfaces [11]:

̇ρo+Viiρo+iVi+voiρo=iDoiρo+σ.E7

This equation can also be put into a vector form:

̇ρo+Vρo+V+voρo=Doρo+σ.E8

### 2.1 Application to oxygen transport

We consider a biological application of the observable’s mass transport equation to microfluidic arterial oxygen transport. In this case, o=O2. For this, we are required to make a few assumptions:

1. (A1) We first tentatively assume that the arteriole’s surface is stationary and not. This would necessarily imply ViZi=0.

2. (A2) We also consider steady-state solutions, by assuming that the density is not dependent on time. This means that tρO2=0.

3. (A3) In addition, we restrict our studies to microfluidic environments which are in equilibrium. This would imply that the net local density change is zero, or σ=0.

4. (A4) We also assume that the diffusion constant is a constant.

Using the above relations, we reduce the conservation relation to

ivO2iρO2=DO2iiρO2.E9

Both of the operators can be expanded using the Voss-Weyl formula [11] and restated in terms of partial derivatives with respect to the spatial coordinates, Zi, and the spatial metric tensor, Zij:

1ZjkZiZjkvO2iρo=1ZjkDO2ZiZjkZiZρO2.E10

We assume for a moment that the coordinate system chosen is some general axial coordinate system consisting of two arbitrary coordinates, Z1Z2, and a third coordinate corresponding to the standard z coordinate found in cylindrical and Euclidean coordinate systems.

This forms a three-dimensional coordinate system of Z1Z2z. We first assume that the velocity of the observable is only along the z coordinate. This means that the term on the left is greatly simplified:

1ZjkzZjkvO2zρO2=1ZjkDO2ZiZjkZiZρO2.

We then assume that the velocity of the observable does not depend on the z-coordinate. This means that the particle moving along a streamline parallel to the length of the tube will not accelerate. This will produce the equation:

1ZjkvO2zzZjkρO2=1ZjkDO2ZiZjkZiZρO2.

From here, if we assume the coordinates Z1Z2=xy, then we can simplify greatly to obtain the final form:

vO2zρO2z=DO2ZiZiρO2Z.E11

In addition, we assume that ρO2 does not depend on x. This produces the final equation of

vO2zρO2z=DO22ρO2y2.E12

In addition, for all the above equations, since by Boyle’s law, at a constant temperature, all instances of oxygen density can be equivalently replaced by oxygen pressure.

## 3. Governing equation for oxygen transport

As can be extrapolated from above, the general form of the nonlinear PDE for consideration defining oxygen transport in core region with Poiseuille hemodynamic flow is given by Eq. (13). The boundary conditions are shown in Section 10.1 of the Appendix, the velocity profile is shown in Section 10.2, and Figure 1 shows the geometry of the problem:

vp1HT+vRBCHTKRBCKp1+HbTKRBCdSO2dPO2PO2z=Dp2PO2E13

There is a core region of blood flow with RBCs and plasma and a cell-free region with only plasma flowing. In the plasma region near the wall, the governing equation is as in Eq. (13), without the second term in the square brackets. In general the geometry of the problem can be either a tube or a channel, and the Laplacian is generalized in [8]. In the present work, we confine the problem to a channel flow. The blood plasma velocity is vp, and vRBC is the velocity of RBCs together with plasma in the cell-rich region. The velocity of the RBCs is lower due to the slip between plasma and RBCs. The distribution of RBCs is such that the hematocrit is higher at the center of the channel and lower near the wall. The term dSO2dPO2 is the slope of the oxyhemoglobin dissociation curve and is a highly nonlinear function of oxygen tension PO2 [1, 2, 3]. The dissociation curve is approximated by the Hill equation [4, 6], where an empirical constant N is used and a constant P50 appears which is the oxygen tension that yields 50% oxygen saturation. HbT is the total heme concentration, Dp is the given oxygen diffusion coefficient in plasma and has units of μm2/s and KRBC, Kp are the solubilities of O2 in the RBCs and plasma, respectively, and have units of M/mmHg. The values of the respective parameters are taken from [8].

In connection with [8], we let

ΣPO2=const+coeffdSO2dPO2.

This allows for a transformation of Eq. (13) into Eq. (14) where “const” and “coeff” are constants derived in [8].

We thus have a PDE of the form:

ΣPO2PO2z=1vp2PO2.E14

Choosing a separation form of PO2

PO2=PO2Z˜z,E15

where Z˜=Z˜1Z˜2 indicates a semi-general coordinate system; we assume the following form of the solution:

ΣPO2Z˜z=PZ˜L1z+L2z.E16

Let Z˜ij be the metric tensor of the coordinate system composed of Z˜1Z˜2. As derived in [8], we express the Laplacian of an arbitrary function, ψZ˜, in terms of the metric tensor in curvilinear coordinates, i.e.,

iiψZ˜=2ψZ˜=1Z˜jkZ˜iZ˜jkZ˜iψZ˜.E17

Applying this definition to 2Σ1PL1+L2 similar to [8], we show that

PL1+L2zΣ1PL1+L2=1vpL1dΣ1dΣ2P+L1Z˜ijiPiPd2Σ1dΣ2.

Rearranging in terms of dΣ1dΣ, we obtain

PL1+L2PdL1dz+dL2dzvpL12PdΣ1dΣ=L1Z˜ijiPiPd2Σ1dΣ2.

Isolating Σ terms and allowing a zero separation constant, we obtain

PL1+L2PdL1dz+dL2dzvpL12P=0.E18

where L10 and P must obey the metric condition

Z˜ijiPiP0.E19

Also, Σ1 must obey the conditions

dΣ1dΣ0d2Σ1dΣ2=0.E20

It has been shown in [8] that for L2z=0, Eq. (13) can be reduced to

2P+S2vpP2=0dL1dz+S2=0,E21

where S2 is a separation constant, vpZ˜=cdr2 is Poiseuille flow, and

PO2=A1PL1+A2,E22

where A1,A2 are arbitrary constants.

In the present work, we consider the general case where the forms are chosen for L1 and L2:

L1=S1Cz+13mi2L2=C1zL1,

where C, C1, and mi are constants. C1 is a free constant, C is to be determined using boundary conditions, and mi is a fixed known constant. The reason for this choice of functions L1z and L2z will be made apparent in Sections 4 and 5. The constant mi is chosen as in Figure 1, and S1 is a free constant.

## 4. Transformation of associated equation

The equation to solve is a PDE related to Eq. (18) and condition Eq. (20), for L1 and L2 defined above:

Prr=S1cdr2{2CCz+1/3miP2+{2zC1Cz+1/3mi+ C1Cz+1/3mi2+2zC1Cz+1/3miC}P+ zC1C1Cz+1/3mi2+2zC1Cz+1/3miC.E23

Let

ξ=ξr=cdr2,E24

The transformed equation becomes

2dPξ+4dcξ2Pξ2=S1ξ{2CMP2+2zC1M+C1M2+2zC1MCP+ zC1C1M2+2zC1MC,}E25

where M=Cz+13mi, z=M13miC, and C1=S11/n=ϵ, small, for n= 5, 4, 3, 2.

## 5. General form of nonlinear equation

The following form of the nonlinear nonhomogeneous PDE, Eq. (25), is considered:

2Yξzξ2+F1ξYξzξ+F2ξzY2ξzϵ/2F3ξb0zYξz+1/2Gξzϵ2=0,E26

where

F1ξ=2c2ξ1,E27
F2ξz=1/4S12CMξdcξ,E28

and Gξzϵ2 involve a small parameter by choice of L1 and L2 and can be made arbitrarily small, whereas F3ξb0z can be made large due to choice of b0z.

In light of work in [7], upon substitution of F1 and F2 for λ as defined in Eq. (33):

λξz=11/4S12CMξdcξ512c2ξ5,E29

and

Yξz=λξzWZξzb0z+ϵMcdt2×z/39Cz+mit2×F31ξb0z,E30

where Z is defined such that

Z=ϕξzξ,E31

where F31=1/F3 is defined by Eq. (35), ϕξz is defined by Eq. (34), and WZ is defined by Eq. (32). It follows that substitution of Eq. (30) for Y into Eq. (26) will cancel part of the coefficient of the Y term leaving the first term on the right side of Eq. (30). (The second term on the right side of Eq. (30) is the simplification of the last term on the right-hand side of Eq. (25) except the F31 term). Next, ε2 terms will cancel in Eq. (26) leaving a form of the equation which is homogeneous similar to that in [7] (where we have assumed that supzb0z is large for z far down the channel), i.e.:

2Yξzξ2+F1ξYξzξ+F2ξzY2ξz+F3ξzYξz=0.

It can be verified that as S10 and ϵ0, the solution obtained from Eq. (26) is a linear function in ξ and independent of z, and one possible solution is a decreasing function on the height of the channel, which is intuitive as the PO2 should drop at the wall of the channel (see Figures 1 and 2). In an approximating sense, though, the other terms will contribute in Eq. (26) as shown below.

The functions λξz,ϕξz and function b0z chosen to be large in magnitude in the supremum sense for all z values are selected so that the transformed equation, Eq.(26), through Eqs. (29)(31) is written as one of the Painlevé classifications of the second-order differential equations [7].

There are two independent canonical forms for this equation [7], one of which is

d2WZdZ26W2Z+6b02=0.E32

According to Estevez et al. [7], the functions must be of the form

λξz=F21/5ξze25F1ξ,E33

and

ϕξz2=λξzF2ξz6.E34

The functions F1,F2 and F3 should satisfy the compatibility relation [7]:

F3ξb0z=2b0zF2ξzλξz2ξ2λξzλξzF1ξξλξzλξz.E35

The previous equation, Eq.(35), after substitution of Eqs. (27)(29) reduces considerably to the following for F3:

F3ξb0z=50cξξ21×12c7ξ+25b0zξ225S1ξCMdcξ4/52c2ξ4/5.E36

The two forms of the solution for F3 (one being the coefficient of P simplified in Eq. (25) and the other Eq. (36)) are equated to each other and compared:

1/5012c7ξ+25b00zM22/5ξ2S1ξCMdcξ4/52c2ξ4/5S11/5CM4/5×
cξ1ξ2=MS1ξcξz2+3C+mi2MS1ξcξz2+3C2.

For z large mi/2 is dropped from the total expression. The ξ4/5 term is scaled by multiplying by a factor of 4 to obtain ξ approximately (see Figure 3). Here b00=1/82+3C/2S11/5

b0z=b00zMS11/5CM4/5.E37

Since z can be large downstream, even for some small S1, then the term 12c7ξ drops out. This results in equality of the two forms of F3 presented. A final substitution for z with a large number multiplied with itself maintains the equality (i.e., zαz, for α large and positive). The solution is valid for z downstream, and we exclude the interval 0mi in Figure 1 with no inlet PO2 value specified as a boundary condition at z=0 as this would give erroneous results in the upstream region.

The final form of the general solution to Eq. (25), downstream (defined on the interval mimf in Figure 1) using Eq. (30) is

Prz=2CM{Nrϵ×z9Cz+mi/6C×F31ξb0z+2CM1π24sinπ2Krz+Q2+1324/5M1rz104zS11/52+3C/16C4/5M1/52CM11M1rz},E38

where

Krz=0.29cdr2CMS1cdr22CMS1cdr21/5,E39
M1rz=2CMS1cdr2)1/5,E40
Nr=cdr2r2,E41

C and Q are constants and ϵ=S11/n, small, for n = 5, 4, 3, 2.

As an approximation, if we divide Eq. (26) by z, using Eq. (30), for large z downstream in the channel, it can be seen that the following equation emerges:

2ξ2λW+F1ξξλW=F2ξλ2W2F3ξλW.

Now the z part appearing in λ is 1/2CM1/5. Multiplying the equation above by 2CM6/5 will result in the left side of the equation to consist of two 2CM terms, and the right-hand side of the equation to have two 2CM9/5 terms where F3ξ has a 2CM4/5 term from Eq. (36). In Eq. (38) the sinπ/2Kr+Q2M1rz term will oscillate to zero as z approaches minus infinity.

A best line of fit can be made by scaling the term 2CM9/5. Multiplying this by approximately 0.01040 results in a best approximation as shown in Figure 4. This allows us to cancel almost all z dependence in the equation except oscillating term and get a homogeneous “almost” z-independent equation as in [7]. Hence, the PDE of Eq. (26) can be reduced to an ode in ξ of similar form but homogeneous as z approaches minus infinity.

## 6. Boundary and matching conditions at wall and core-plasma interface

We employ the Robin boundary condition Eq. (58) in Appendix and derivative matching conditions at the interface of plasma and RBC core regions, respectively, at z=mi and z=mf, shown in Figure 1, to determine constants Q and C as well as two additional constants for linear solution in plasma layer. The solution to the linear part of Eq. (13) defined in the plasma layer, i.e.,

vpPO2plasmaz=Dp2PO2plasma,E42

is

PO2plasma=CP1CylinderU1/2cS1Dp1dr4S1dDp4+
CP2CylinderV1/2cS1Dp1dr4S1dDp4,

where CylinderU and CylinderV are parabolic cylinder functions and CP1 and CP2 are constants to be determined. S1 is a separation constant for Eq. (42). The following boundary matching condition is utilized at the interface of the plasma layer (1 micron in height) and RBC core region (49 microns in height):

PO2plasmar=PO2corer,E43

at z=mi and z=mf, respectively (see Figure 1). Four equations in four unknowns were solved in Maple 18, where two constants are from each of the two solutions in plasma and core regions, respectively. The total of eight constants was determined and is shown in Table 1 for low and high hematocrit. The value of Dp is obtained from Table 1 in [6] and c=1250 and d=0.5. For high hematocrit as shown in the Appendix, the solution incorporates different values of c and d. Also we let S1=0.911×108 in the core region and S1=0.211×103 the in plasma layer. Here there are two different separation constants for two different regions.

ConstantsCQCP1CP2
High hematocrit = 0.45−0.10096871926.261616672 ×1071.1113336191.542724257
Low hematocrit t = 0.15−0.098247473182063.8319661.1113336191.542724258

### Table 1.

System of constants for associated system of four unknowns solved in maple.

The fourth equation used was that the change in flux at the wall at the edge of the plasma layer far downstream z<<mf is zero. In addition the no-flux condition shown in Figure 1 (i.e., impermeable membrane) at r=0 was satisfied exactly for all z for the core region solution. This is the fifth boundary condition. There is no inlet PO2 specified at z=0. The oxygen tension in core is PO2core=33.07440049+3.638040582107P which is in the form of Eq. (22) and gives a PO2 of approximately 150 mmHg at r=0, z=mi=7000 microns.

## 7. Weierstrass elliptic function

The Weierstrass P function is defined as

WZ=1Z2+w1Zw21w2.E44

As in [9], we consider the following expression:

η=expuπi/ω,E45

and derive a function of η which behaves like the Weierstrass P function at η=0.

The development of η in the neighborhood of u=0 is

η=1+uπiω+12!uπiω2+,E46

or

η1=uπiω1+12!uπiω+13!uπiω2]E47

Observe that the function

η12=u2π2ω21+uπiω+712uπiω2+,E48

is zero of the second order at all the points u=2μω(μ=0,±1,±2,±3,..).

Consider a function Jη such that Jη0 for η=1, the function

Jηη12,E49

is infinite of the second order for all values u=0 and u=2μω.

This behavior at these points is the same as the Weierstrass P function.

We write Jη=a++cη2 where a,b,c are constants.

Since η2=e2uπiω, we have

Jηη12=a+b1+uπiω+12!uπiω2++c1+2uπiω+12!2uπiω2u2π2ω21+uπiω+13uπiω2+ω2π2a+b1+uπiω+12uπiω2++c1+2uπiω+u21uπiω+u2.E50

As in [9] the Weierstrass P function is shown to have the following degenerate form which is used in Section 5:

Jηη12=π/2ω2sin2/2ω1/3.E51

## 8. Results and discussion

The present work shows that near the inlet of the permeable membrane (see Figure 1), there is a significant drop at the wall of the channel in PO2 as compared to downstream values as shown in Figures 2 and 5. It is worthy to note that the structure of the solution obtained in terms of degenerate Weierstrass P function forms a near constant solution across the height of a micro-fluidic channel downstream at end of the membrane region. This is seen in both Figures 2 and 5, and in Figure 5, the contours flatten out as the flow of blood proceeds far downstream. The release of ATP has been shown to be caused by a change in oxygen saturation [6]. It is concluded that there is a significant decrease in oxygen tension to the right of the permeable membrane downstream. It is in this region that there is a significant concentration of ATP released. It is important to mention that a rapid decrease in oxygen saturation is one means to produce ATP; it is also accomplished by means of applying shear stress on the RBC as shown in [12].

## 9. Conclusion

A well-known governing nonlinear PDE used to model oxygen transport was formulated in a recent paper in a generalized coordinate system where the Laplacian is expressed in metric tensor form. A reduction of the PDE to simpler problem subject to specific integrability conditions was shown there. A reduced almost linear ode was derived, and in the present paper, a solution has been obtained using a well-known factorization method for the second-order ode where a compatibility equation has been used in equating it to a specific form of the original differential equation. Approximate oxygen tension profiles have been determined downstream in a micro-channel in the vicinity of a permeable membrane with an oxygen supply on the other side of the membrane. Although it is expected that ATP will be released as blood flows past the permeable membrane downstream, it has been shown mathematically that this is the case and increases in hematocrit produce more ATP. Future work remains to apply tensor equations for a moving arterial surface as generalized at the start of the present work.

## Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

## A.1 Cell free

For cell-free region as shown in Figure 1, we have that

vpPO2z=Dp2PO2,E52

where vp is the velocity of plasma, given at the end of the Appendix.

## A.2 Core region

In this case there is a central core region where oxygen dissociates to form HbO2.

Therefore, the velocity in core region, i.e., vO2=fvpqidSO2dPO2.

The model consists of the following partial differential equation:

vp1HT+vRBCHTKRBCKp1+HbTKRBCdSO2dPO2PO2z=Dp2PO2.E53

## A.3 Boundary conditions

The boundary conditions are defined as follows:

PO2rmi=PO2i,r0h,E54
1. Pre-membrane

PO2rr=h=0,z0mi.E55

2. Bottom of membrane region

PO2rr=0=0,z0L.E56

3. Post-membrane

PO2rr=h=0,zmfL.E57

The only region not covered by the above three flux equations is the region of PO2 occupying zmimf at r=h the membrane. This region is governed by a Robin condition specified in [8].

The boundary condition is

PO2rr=h=DmKmDpKpPO2oPO2hzτ,zmimf,E58

where PO2o is the PO2 level on the other side of the membrane.

All of these conditions act on the entire system, and constants appearing are found in Table 1 in [6].

## A.4 Velocity profile function

The following velocity profile is used in channel:

vpx̂=31301064wh3+μp/μcyi3h2x̂2,yix̂hh2yi2+μp/μcyi2x̂20x̂yiE59

Relative apparent viscosity, μp/μc1, for low discharge hematocrit. The previous results (Figures 2 and 5) were based on this data. Figure 6 is based on a discharge hematocrit of approximately 40% higher with relative apparent viscosity equal to 1.7. It is apparent from the graph that in comparison to Figure 2 with lower discharge hematocrit that there is an increase in the drop of PO2 profiles downstream with a resulting higher production of ATP.

## References

1. 1. Nair PK, Huang NS, Olson JS. A simple model for prediction of oxygen transport rates by flowing blood in large capillaries. Microvascular Research. 1990;39:203-211
2. 2. Nair PK, Hellums JD, Olson JS. Prediction of oxygen transport rates in blood flowing in large capillaries. Microvascular Research. 1989;36:269-285
3. 3. Nair PK. Simulation of Oxygen Transport in Capillaries. Rice University; 1988
4. 4. Moschandreou TE, Ellis CG, Goldman D. Influence of tissue metabolism and capillary oxygen supply on arteriolar oxygen transport: A computational model. Mathematical Biosciences. 2011;232(1):1-10
5. 5. Ng C-O. Dispersion in steady and oscillatory flows through a tube with reversible and irreversible wall reactions. Proceedings of the Royal Society A. 2006;462:481-515
6. 6. Sove RJ, Ghonaim N, Goldman D, Ellis CG. A computational model of a microfluidic device to measure the dynamics of oxygen-dependent ATP release from erythrocytes. PLoS One. 2013;8(11):1-9
7. 7. Estevez PG, Kuru S, Negro J, Nieto LM. Factorization of a class of almost linear second-order differential equations. Journal of Physics A: Mathematical and Theoretical. 2007;40:9819-9824
8. 8. Afas KC, Moschandreou TE. Analytic multiplicative separation and existence investigation of non-linear oxygen transport with Poiseuille hemodynamic flow in a semi-generalized co-ordinate system. International Journal of Differential Equations and Applications. 2015;14:281-311
9. 9. Hancock H. Lectures on the Theory of Elliptic Functions. Dover; 2004
10. 10. Brennan MD, Rexius-Hall ML, Elgass LJ, Eddington DT. Oxygen control with microfluidics. Lab on a Chip. 2014;22:1-14
11. 11. Grinfeld P. Introduction to Tensor Analysis and the Calculus of Moving Surfaces. Springer; 2010
12. 12. Wan J, Ristenpart WD, Stone HA. Dynamics of shear-induced ATP release from red blood cells. PNAS;105(43):16432-16437

Written By

Terry E. Moschandreou and Keith C. Afas

Reviewed: November 20th, 2019 Published: December 26th, 2019