Following red blood cells in a pulmonary capillary
Abstract
The red blood cells or erythrocytes are biconcave shaped cells and consist mostly in a membrane delimiting a cytosol with a high concentration in hemoglobin. This membrane is highly deformable and allows the cells to go through narrow passages like the capillaries which diameters can be much smaller than red blood cells one. They carry oxygen thanks to hemoglobin, a complex molecule that have very high affinity for oxygen. The capacity of erythrocytes to load and unload oxygen is thus a determinant factor in their efficacy. In this paper, we will focus on the pulmonary capillary where red blood cells capture oxygen. We propose a camera method in order to numerically study the behavior of the red blood cell along a whole capillary. Our goal is to understand how erythrocytes geometrical changes along the capillary can affect its capacity to capture oxygen.
The first part of this document presents the model chosen for the red blood cells along with the numerical method used to determine and follow their shapes along the capillary. The membrane of the red blood cell is complex and has been modelled by an hyperelastic approach coming from [16]. This camera method is then validated and confronted with a standard ALE method. Some geometrical properties of the red blood cells observed in our simulations are then studied and discussed.
The second part of this paper deals with the modeling of oxygen and hemoglobin chemistry in the geometries obtained in the first part. We have implemented a full complex hemoglobin behavior with allosteric states inspired from [4].
Les globules rouges ou érythrocytes sont des cellules biconcaves et sont formées d’une membrane contenant le cytosol, fluide ayant une grande concentration en hémoglobine. Cette membrane est très déformable et permet à ces cellules de traverser des passages étroits comme les capillaires dont les diamètres sont souvent plus petits que ne l’est le diamètre des globules rouges. Les érythrocytes transportent l’oxygène dans l’organisme grâce à la proteine d’hémoglobine qui dispose d’une grande affinité pour l’oxygène. La capacité des globules rouges à capturer et libérer l’oxygène est ainsi un facteur déterminant de leur efficacité. Dans cet article, nous nous restreindrons à l’étude d’un capillaire pulmonaire où les globules rouges capturent l’oxygène. Nous développons une méthode de caméra de façon à étudier numériquement le comportement de globules rouges le long d’un capillaire. Notre but est de mieux comprendre comment les éythrocytes se déforment au long du capillaire et comment cela affecte leur capacité à capturer l’oxygène.
La première partie de ce document présente le modèle que nous avons choisi pour représenter le globule rouge, ainsi que la méthode numérique utilisée pour déterminer et suivre la forme du globule le long du capillaire. La membrane du globule rouge est complexe et a été modélisée par une approche hyperélastique issue de [16]. Cette méthode de caméra est ensuite validée et confrontée à une méthode ALE standard. Enfin, quelques propriétés géométriques du globule rouge résultant de nos simulations sont étudiées.
La deuxième partie de ce papier traite de la modélisation de la chimie entre l’hémoglobine et l’oxygène dans des géométries issues de la première partie. Un modèle complet du comportement de l’hémoglobine a été implémenté, en particulier faisant intervenir deux états allostérique des molécules d’hémoglobine. Ce modèle s’inspire de [4].
Introduction
Human respiration and more generally mammal respiration, is a complicated efficient system that uses oxygen amongst other products in order to provide energy to the body. One of the function of the respiration is thus to bring oxygen from air into the mitochondria where the major energy source of the body, the ATP, is synthesized. Oxygen travels a complex exchange system in different media. Ambient air is first brought on an air/blood exchange surface (the aveoli) through the lung, where oxygen is extracted from air. The fractal like structure of lungs (dichotomous tree of twentythree generations) permits the exchange surface to be large enough () to fit energy needs of the body. The particular structure and functioning of the lung are very important in respiration and have been intensively studied, see for example [9, 15] but mostly independently of blood behavior. The next step of respiration for oxygen is to cross the exchange surface, to reach the red blood cells in the capillaries (smaller blood vessels) and to react with hemoglobin in order to be carried away from the lung towards the cells and mitochondria through an other branched complex structure, the vascular network, which also has been studied intensively, see for example [10].
In this paper, we are interested in the modeling of the step of oxygen transport that takes place in the pulmonary capillaries. There, oxygen crosses the membrane between the alveola and the capillary and dissolves into plasma, where oxygen is transported by a convection/diffusion process towards the red blood cells. Next, oxygen crosses the membrane of the red blood cell and diffuses inside it (in the cytosol), then meeting an hemoglobin molecule, a chemical reaction can occur and oxygen is trapped by the hemoglobin [23].
Such a modeling can lead to a better understanding of the phenomena involved in oxygen exchanges at capillary level. In particular, some pathologies can affect these exchanges by modifying critical parameters. For example, the rigidity of the red blood cells membrane can be modified (paludism, diabetes) and lead to less efficient shapes in terms of capacity to cross capillaries or capacity to catch oxygen. Similarly, the chemical properties of hemoglobin are very sensitive to blood parameters (as pH) and global consequences on oxygen transport in blood could be determined.
The convection/diffusion process involved in oxygen transport in the capillaries is highly dependant on the shape of the structures (capillary and red blood cells) but also on the velocity field of the plasma. Hence the determination of these are required in order to correctly model oxygen displacements. The capillaries diameters are often much smaller than the red blood cells’ ones, but thanks to their highly deformable membrane the red blood cells shapes can change under external constraints and they can cross the smaller vessels (parachute shape).
In order to deal with these different points, a numerical camera method is used to follow the red blood cells deformation along a straight capillary. This method is based on the standard ALE method for fluid structure interaction but is written in a moving frame. The benefits of this approach is that it is not limited by mesh skewing and the whole capillary can be covered. The mechanics of the membrane of the red blood cells is modeled by a Yeoh hyperelastic model from [16] and validated. Some properties of red blood cells (velocities, shapes in the capillaries) are studied along with apparent hydrodynamical resistances of capillaries crossed by red blood cells.
The second step of this paper is to add to the previous modeling the behavior of oxygen in the plasma and in the cytosol, and the behavior of hemoglobin in the cytosol. Two physical processes are involved: displacement (convection and diffusion for oxygen, diffusion only for hemoglobin) and chemistry (oxygen and hemoglobin reactions). A realistic model of hemoglobin chemistry is used [4] and implemented numerically. Some preliminary results are presented, mostly on the consequences of capillaries diameter on the red blood cells capacity to capture oxygen.
Note that all simulations have been performed using Comsol Multiphysics 3.3a and that most of our calculations are 2Daxisymmetric.
1 Fluid structure interaction: modeling
1.1 Geometries
1.1.1 Red blood cells
The red blood cell can be considered as vesicle with a membrane enclosing a fluid, the cytosol. The shape of the red blood cell is particular and known to minimize the bending energy of the membrane [2] amongst other properties [22], see figure 1. In the following work, the geometrical parameters of the red blood cells have been chosen accordingly to Weibel’s data [23] for human erythrocytes.
The cytosol can be considered as a viscous fluid and contains a high concentration of hemoglobin. It will be considered as a Newtonian viscous fluid: its density corresponds to that of water while its viscosity is eight times water viscosity [14].
The membrane of the red blood cells is a superposition of a bilipidic layer and of an actin mesh. Its thickness is of the order of one or two hundreds of nanometers [23].
The mechanical behavior of the membrane is complex and not yet well understood. However some interesting models have been developed [1, 16, 18]. More precisely, in [16], a relatively simple hyperelastic model for the membrane has been applied successfully to recreate numerically optical tweezers experiments. This hyperelastic model has been developed by Yeoh in [25] and uses the following strain energy:
(1) 
where is the initial shear modulus of the membrane and a constant estimated to in order to fit optical tweezers experiments, see [16]. Note that in the literature, data are often given in term of inplane shear modulus: where is the initial thickness of the membrane. In the following, has been chosen to , see 1.1.2.
1.1.2 Reproducing Mills et al experiments
In order to validate the membrane model, a numerical reproduction of Mills et al experiments [16] has been performed. These experiments consist in attaching a metallic bead at each extremity of the red blood cell and to use magnetic sinks in order to control the positions of the beads (optical tweezers). This method gives access to the force applied to the beads and thus it is possible to deduce some mechanical properties of the membrane if measuring the deformations against the force applied.
In order to build a numerical equivalent, we use our three dimensional red blood cell. Two cubic rigid “beads” have been attached on its membrane. One of these beads will be fixed while the other will be submitted to a force. The inside of the red blood cell is considered as a Stokes fluid with previously given characteristics. Note that because the fluid is assumed incompressible, the red blood cell keeps a constant volume. We use an ALE fluidstructure interaction model. Figure 2 shows the initial geometry and a typical deformed globule under stress.
We have measured lengths, widths and thicknesses of the red blood cell for a range of forces going from to with steps of . We have found a better collusion between our results and Mills et al results for an inplane shear modulus of and with the parameter equal to . The corresponding results of our numerical simulations are shown on figure 3. In the following, we will keep these parameters constant to these validated values.
1.1.3 Capillary
We will assume at first approximation that a capillary is topologically equivalent to a cylinder whose axis can be considered as a straight line. In the camera method developed below, this will simplify the movement of the camera, limiting it to a translation along the capillary axis. The walls of the capillary will be assumed unmoving and rigid. The shape of the capillary will be chosen according to measures done in [21]. The typical length of such a vessel is of the order of one millimeter and its radius can range from to . Figure 4 shows a geometrical model of capillary that will be used in the next sections.
The plasma will be considered as a viscous Newtonian fluid with same density than water, but with twice its viscosity. Because capillaries are far from the heart, they do not feel heart beating and thus we can consider that the plasma circulates with a constant flow rate along the capillary. Moreover the velocities are of the order of the millimeter per second and the Reynolds number can be estimated to . Thus Stokes equations are a good approximation of plasma behavior in the capillary.
1.2 Equations
Because red blood cells are travelling through the whole capillary, the coordinates map induced by their displacement in a standard ALE method will be highly skewed. Direct numerical simulations using finite elements will thus be very difficult because of highly deformed elements. This could be bypassed by regular remeshings but this method is relatively costly and relatively difficult to set up.
Hence, considering the typical shapes of the studied capillary (that is assumed to be topologically equivalent to a cylinder and to keep its axis along a straight line) and considering the fact that our objects of interest are the red blood cells, a camera method seems advantageous. This method consists in focusing only on the portion of the capillary where the red blood cells are and to displace focus together with red blood cells.
Numerically speaking, this consists in meshing a small portion of the capillary, in making the mesh move along it at red blood cells mean velocity and also in making the mesh fit capillary local shape. Note that because we do not need to mesh the complete capillary, this method will also benefit of a reduced number of elements.
The drawbacks of this method are located on boundaries:

the coordinates map induced by red blood cells deformations must also include geometrical changes of the capillary walls along its length

the boundary conditions for the plasma on the nonwalled boundaries are not obvious and need to be studied more closely
Note that all bearings involved in the following will be built thanks to the Winslow smoothing method [24], which robustness against mesh elements inversion is better than Laplace method.
1.2.1 Standard method
First we will write the standard method, namely the standard fluid structure equations in ALE formulation on the full capillary. The capillary domain will be called and its straight axis will correspond to one of the reference frame axis. In the following, this particular direction will be recognized amongst the other by the subscript . The reference coordinate system will be denoted , and the basis vectors will be noted . Thus the cylinder axis is directed by and the position on this axis referenced by .
Let be the membrane of the red blood cells and the reference domain for the fluid. Note that corresponds not only to the outside of the red blood cells (plasma), but also to their insides (the cytosol), these two fluids have the same density (similar to water) but a different viscosity. Let () be the capillary walls as shown on figure (5) such that .
The standard fluid structure equations in ALE formulation are:
(2) 
with boundary conditions:
(3) 
and , where is a bearing of on whole with boundary constraints on and on . The constraints equality is done on the reference frame and the notation means that the constraints tensor of the fluid defined in the deformed frame is pulled back in the reference frame. The boundary condition on is the displacement due to the deformation of the membranes of the red blood cells, this deformation comes from equations (2). The boundary condition on states that the capillary walls are unmovable but allows the coordinates of the deformed frame to slide along the walls of the capillary. The coordinate in the deformed frame will be noted and .
1.2.2 the camera method
We will now restrict the domain for the camera method: we will now consider a coordinates system directed by and assume our reference frame is a cylinder which axis is along . Its length along is ( is smaller than capillary length). This cylinder will fit the local capillary geometry while it follows the red blood cells along the capillary. The deformed cylinder at time will be called the camera and will be noted . The domain defined by the red blood cells membranes is called .
The camera will follow the red blood cells in such a way their mean position is always at the center of the view of the camera. Hence the position along the capillary of the camera center can be calculated:
(4) 
The displacements are still calculated in the reference frame on the membrane domain and will not differ from the standard ALE method. In order to rebuild the geometry of the capillary containing the deformed globules, it is necessary to first integrate the deformed globules into the nondeformed capillary and then to make the cylinder wall fit the capillary without changing the shapes of the globules. This will be done using two successive bearings and :
(5) 
Thanks to the definition of , the transformation preserves the mean position of the red blood cells to on the axis.
Stokes equations (velocity will be called and pressure ) will be calculated in coordinates in order to be solved in the correct geometry. The transformation corresponds to the movement of the camera and will thus induce an additive transport term in the Stokes equations: . However, because of the small Reynolds numbers involved and thanks to the fact that , this supplementary term is negligible the same way is. Moreover, the transformation does not modify Stokes equations.
Finally, the succession of these two bearings can be reduced to a single bearing . The boundary conditions for are:
(6) 
The boundary conditions for fluid equations will be the same as for the standard ALE method concerning the capillary wall and the globules wall (non slip conditions). However, they should be studied for what concerns the inlet () and the outlet () of the camera domain. In the standard ALE method, a Dirichlet condition on velocity was imposed at the outlet and a zero pressure condition was introduced at the inlet . The Dirichlet condition can be easily used at camera outlet considering the very low Reynolds numbers: a Poiseuille profile should be a good approximation of the velocity at camera outlet, all the most if the border is not too close to the red blood cells. At inlet () however, the correct pressure condition should be with such that and such that . Again, thanks to the low Reynolds numbers, we can assume that pressure is constant on the section and impose to fit the mean value on this section namely . However this value is only time dependant and because Stokes equations only depends on , making the variable change will not change the equations and will allow us to impose on the camera inlet.
Finally, the fluid constraints on red blood cells membranes is modified only by the relative pressure change: because this change is global it does not affect the mechanics equations.
The errors done relatively to the standard approach are due to boundary conditions simplifications and to the suppression of the new transport term. We will do numerical comparisons in order to show the quasi equivalence between both methods and show that, in the contrary of the standard approach, the camera method can follow the red blood cells all along the whole capillary.
2 Fluid structure interaction: numerical results
2.1 Comparison between the standard method and the camera method (2D axisymmetric)
In order to validate the camera method, a comparison has been made between the standard ALE method and its pendant in camera method. In order to limit the calculation time, a 2D axisymmetric case has been used. In this example, a red blood cell has been put inside a plasma flow through an axisymmetric bottleneck. At first the radius of the capillary is larger () than that of the globule () but it is smaller in the constriction (), see figure 7.
Because of its design, the camera method induces a far smaller distortion of the mesh and can follow easily the red blood cell along the whole geometry. The mesh quality for the camera method is only dependant on the radius of the capillary and not on the position of the red blood cell. On the contrary, mesh quality is decreasing linearly in the standard method, and the numerical simulation stops before the end of the bottleneck (at time about from entry) because of highly skewed elements. Mean mesh qualities versus globule positions are represented on figure 8.
Some parameters have been checked in order to compare the differences between both methods. We have measured the displacement of the red blood cell and also its aspect ratio along time in both simulations. The data show a very good agreement between both methods on the range of convergence of the standard method. There is however a slight difference for the red blood cell aspect ratio in the last times of convergence of the standard method that is a consequence of some of its elements being highly skewed: this lead to a bad convergence of the problem solution.
2.2 Going along a capillary
Being able to follow red blood cells along a capillary gives access to a large quantity of data that needs to be studied in lots of different contexts. We will present here some of these interesting results but keeping in the wide context: “how do it feel being a red blood cell in a capillary ?”.
In this section, and as before, we will limit us to a 2D axisymmetric problem in order to limit calculations times.
We will do our study on an infinite capillary along axis (measures will be done in meters). The vessel diameter is equal to on and is variable on the range . The variations chosen are close to the one observed on an invivo rat capillary in [13] and are shown on figure 10. We will call the portion of capillary corresponding to the region of diameters variation, namely with .
Depending on the case, we will put one, three or five red blood cells in the capillary. The initial distance between them has been put to in order to avoid collisions. This value is higher than the observations that gives a number closer to (data calculated from a red blood cells volume fraction of in blood, see [19, 23]). So that they have a stabilized shape at the entry of the portion of capillary , the initial mean position of the globules has been put to . Plasma flow through a capillary section of diameter has been chosen to correspond to a mean velocity of , in the simulations the flow is maintained constant whatever the position of the camera, in particular it implies velocity increase when capillary radius decreases.
2.2.1 What shape?
Under plasma flow stress, the red blood cells are moving forward the capillary, and because of the non homogeneous shape of the flow on a capillary section, the globules shape changes to the typical parachute shape, pointing in the direction of the flow. Moreover the globules go through skewer shapes when the radius of the capillary decreases. In order to estimate the distortion of red blood cells, we will use the aspect ratio as in [13]. It is defined by the ratio of cells length along the capillary axis over their diameter, see figure 11. Under no stress, this aspect ratio is .
Figure 11 shows the different curves of aspect ratio depending on the number of red blood cells involved in the calculation. Red blood cells aspect ratio follows radius changes of capillary. This illustrate how able to deform the red blood cells are: starting from an aspect ratio of , they can reach almost in the narrowest parts of the capillary. These numbers are fully coherent with the measures made in [13], where the authors have observed similarly the red blood cells aspect ratios and have more particularly studied these changes against capillary radius variation. On figure 12, typical shapes observed in our simulations have been drawn, in particular we can observe the classical parachute shape.
The presence of more than one red blood cells plays a small but visible role in aspect ratio determination for our choice of parameters (distance of between each cell). This indicates that closer globules, as in real blood, should influence themselves, all the more if there are collisions between them.
2.2.2 Which velocity?
Because the skewed red blood cell does not fill the whole section of the capillary, it is submitted only to the larger part of plasma velocities which stands mostly in the middle of the capillary (Poiseuille profile). Thus, the mean velocity of the cells is most of the time larger to the mean velocity of the plasma. This characteristic of blood is well known [7]. Similarly, it is interesting to observe that the gap between the cells and the capillary wall decreases with capillary radius and can indeed modify this phenomena. The velocities of plasma and the velocities of cells in different cases (one, three or five cells in the capillary) are shown on figure 13. Note that, for the left figure, increasing the number of cells in the capillary will induce an homogenisation of the mean velocity because each cell sees a different radius than the others. On the right part of the figure the velocities of the middle red blood cell is shown, the effects of the presence of other cells is small but exists. Hence, decreasing the initial distance between cells should induce a higher interaction, but also create collisions which are not, for now, managed by our code.
2.2.3 Consequences on the hydrodynamical resistance of the capillary?
By crossing the capillary, the red blood cells interfere with vessel’s hydrodynamical resistance which cannot be predicted through fluid models, even nonnewtonian ones. This is due to the deformability of the red blood cells along with the similar sizes of red blood cells and capillary diameters. Thus it is necessary to study the full fluid structure interaction problem in order to obtain reasonable high scales estimation of capillaries hydrodynamical resistances. In order to apply our simulations at this problem, we have obtained preliminary results on the hydrodynamical resistance. We have measured the resistance of the portion of capillary viewed by the camera against the number of cells (zero, one, three or five). To estimate resistance we have calculated the ratio from pressure drop between the inlet and the outlet over the flow in a capillary section.
The results are drawn on figure 14 and show that, as expected, increasing the number of cells increases the resistance. Note that at relative high diameters, the resistance is almost identical whatever the number of the globules. On the contrary, the differences increases when globules shape is changing, actually some of the fluid energy is exchanged with the membrane of the red blood cells and higher pressure drop is thus needed to ensure constant flow.
3 Introducing oxygen and hemoglobin: modeling (preliminary)
Red blood cells are designed to catch oxygen in the lung capillaries, to carry it in the whole body and to deliver it to the cells in the systemic capillaries. The way oxygen is captured and delivered is obviously critical. Actually, the time spent in the capillaries must be tuned in order to maintain the balance between oxygen delivery and oxygen needs. The chemistry of hemoglobin plays an important role in oxygen exchanges, but other actors are also involved, like blood velocity or red blood cells shapes that are determinant for good functioning of the system.
In this section we will use the hemoglobin chemistry model developed in [4] and extend it to axisymmetric one in the geometries obtained from the preceding sections. More precisely and because oxygen diffusion and reaction with hemoglobin can be uncoupled from red blood cells deformation, the results of the preceding sections can be considered as known data in the modeling of oxygen exchanges.
3.1 Modeling of hemoglobin chemistry
Hemoglobin is a complex molecule consisting in four hems, each of them being able to carry one molecule of oxygen. The four hems divides in two and two hems (or chains) that have different bindings properties with oxygen. This asymmetry of the hemoglobin molecule implies the existence of two allosteric states and which will be integrated in the model. It is not clear in which priority these chains react with oxygen, and we will follow one of the mechanisms proposed in [4] which considers that chains are the first to react. With this hypothesis, the rate constants of the different chemical reactions have been determined by multiple experimental measurements which references are available in [4]. The diagram of the reactions involved in the model is shown on figure 15. Note that because these experiments are not made at the temperature of human body, adjustments have been made by multiplication with a unique factor each of the constant rate of bimolecular reactions. Once these adjustments done, we have successfully compared the resulting chemical model with known properties of hemoglobin.
In the following, we will write the concentration of hemoglobin of state ( corresponds to state while corresponds to state) carrying oxygen molecules (noted on figure 15).
An important number in blood studies is hemoglobin saturation, this number corresponds to the percentage of hemoglobin capacity which is used to carry oxygen. For instance, during rest regime hemoglobin saturation is after crossing the lungs while it is before. It is the ratio between the quantity of effectively captured over the maximum quantity of possible to capture. In the model chosen here, it is calculated by:
(7) 
On figure 16, the saturation curve resulting from the model is shown. It has been plotted against oxygen partial pressure in the medium (). The steep part has an important role in oxygen exchange and its position is important, see [23] for more details. The model fits well measurements, for example the partial pressure in oxygen needed to have a saturation (steep part of the curve) is while it is measured to in [23].
3.2 Equations
The path of oxygen from the aveolus into an hemoglobin molecule is described on figure 17. In order to define the equations, we will note the plasma region and the cytosol. We will integrate membranes crossing in the boundary conditions in order to simplify the problem, we will note the boundary between alveolus and capillary and . We will call the plasma velocity which is considered as a data in this part, but which comes from calculations of the preceding sections. We will call the oxygen partial concentration in plasma and the oxygen concentration in the cytosol .
Oxygen equations on are convection/diffusion equations written in the deformed frame for all :
(8) 
where is the diffusion coefficient of oxygen in plasma, the diffusion coefficient of oxygen in the tissues, is the alveolus/capillary membrane thickness, the solubility coefficient of oxygen in plasma and the oxygen partial pressure in the aveolus. The supplementary term corresponds to the fact that the camera is moving and determines when it will reach a region permeable to oxygen. Thus, the set corresponds to the portion of capillary which is permeable to oxygen. Typically, in our simulations : with being the capillary typical length (), see section 4. We assume that no diffusive flow crosses the inlet and the outlet. Finally is the thickness of the red blood cell membrane and we recall that is the concentration of oxygen in the cytosol .
Note that because is close to near the red blood cell, diffusion is the major phenomena involved in the transport of oxygen near the globule. However, close to the wall of capillary, the apparent transport velocity of oxygen is (because ) and this velocity is of the same order than the diffusive velocity ( is the mean distance to be crossed by diffusion). Hence convection and diffusion each play a non negligible role in oxygen repartition near the capillary wall (the Péclet number is ).
We will now write hemoglobin and oxygen equations in the cytosol . Each hemoglobin molecule will diffuse in the cytosol (no transport) and can either:

capture an oxygen molecule except if , at rate .

free an oxygen molecule except if , at rate .

change its configuration, , at rate .
Hence, doing the balance of sources and sinks of each quantities, we can write, in the deformed frame, the following equations for concentration , for all and :
(9) 
is the hemoglobin diffusion coefficient in the cytosol. The boundary condition means that hemoglobin does not cross the red blood cell membrane.
Oxygen diffuses into the cytosol and reacts with the different molecules of hemoglobin. The balance for oxygen in the cytosol can then be written in the deformed frame for all :
(10) 
We assume that the diffusion coefficient of oxygen in the cytosol is smaller than in plasma, because hemoglobin are big molecules. It has been estimated to be of the same order than in the tissue, see [12].
The initial conditions are chosen in order to correspond to poor oxygenated blood with of oxygen partial pressure. Initial hemoglobin state is assumed at equilibrium everywhere in the red blood cell with this oxygen partial pressure.
4 Oxygen and hemoglobin: first numerical results
We will restrict this preliminary study to two similar cases with pulmonary capillaries which geometries are simple bottlenecks as on figure 7, but with different sizes: the structure is long and the constricted part, where oxygen exchanges will take place, is long (typical length of a capillary see [23]). The larger diameter will be in both cases. The smaller diameter will be in the first case and in the second one. These geometries will be crossed by five red blood cells. We assume that oxygen exchanges with alveolus starts at the beginning of the bottleneck.
On figure 18, the mean saturation of oxygen is plotted against the mean position of the red blood cells in the capillary. The first case, with smaller diameter, induces a quicker capture of oxygen by hemoglobin. This is mainly due to the fact that:

the distance between the red blood cells and the capillary wall is smaller in the first case than in the second case, hence oxygen crosses the gap quicker when the diameter is small.

the red blood cell surface “available” for oxygen from capillary wall is larger when the diameter is smaller. Actually, for the larger diameter case, a large part of the surface stays in the middle of the capillary and oxygen needs to cross a larger distance to reach it.
These observations are shown on figure 19, where the middle red blood cell is drawn for both cases. Both images are taken at the same abscissa, and their difference in term of oxygen saturation is important. As expected the red blood cell regions the more accessible to oxygen (closer to the capillary wall) are first saturated even though hemoglobin is diffusing inside the cytosol. Note that the less saturated region is on the axis, in the rear part of the red blood cell where oxygen flow is the smaller. This effect is due to the screening made in this region by the red blood cell shape to the diffusive process. The screening effect of diffusion can thus play an important role in biological systems. This influence has already been shown in the case of the lung, see [8]. Note that the presence of the other four red blood cells also limits oxygen capture speed through the same phenomena.
Conclusion
In this paper, a full modeling chain of red blood cell and oxygen behaviors in a capillary has been built. This gives access to a large panel of data which should be useful in order to better understand oxygen transfer at capillary level. For instance, the fine tuning of some physiological parameters (for example as in the lungs or as blood speed regulation in the vascular system) should be explained by the phenomena involved in capillaries and help to understand the evolution processes that led to such structures. Moreover, pathological behavior (for example paludism or diabetes) could be integrated in this model and could help to understand the consequences of such dysfunctions.
The prospects of this approach is to develop a three dimensional model and to add collisions between red blood cells and with capillary walls. There also exists an important interaction between erythrocytes: they can combine together to form rolls, up to eight cells. These rolls are maintained by chemical forces which can break under high shears. The integration of this behavior is important because of its role in blood circulation. Moreover the restriction to straight capillary could be removed. Finally, mathematical studies of the method should be performed in order to control precisely the errors made during the modeling.
acknowledgement :
The author would like to thank Philippe Dantan from MSC (Paris 7) for useful discussions on blood properties and for his help with Comsol Multiphysics.
References
 [1] Arslan M., Boyce M.C., Constitutive modeling of the finite deformation behavior of membranes possessing a triangulated network microstructure, J. of App. Mech., vol. 73, pp 536543, 2006.
 [2] Canham P.B., The Minimum Energy of Bending as a Possible Explanation of the Biconcave Shape of the Human Red Blood Cell, J. theor. Biol., vol. 26, pp. 6181, 1970.
 [3] Cottet G.H. and Maitre E., A levelset formulation of immersed boundary methods for fluidstructure interaction problems, C. R. Acad. Sci., vol. 338, pp. 581586, 2004.
 [4] Czerlinski G., Levin R. and Ypma T., Hemoglobin/ Systems: Using Shortlived Intermediates for Mechanistic Discrimination, J. Theor. Biol., vol. 199, pp. 2544, 1999.
 [5] Dao M., Lim C.T. and Suresh S., Mechanics of the Human Red Blood Cell Deformed with Optical Tweezers, J. Mech. Phys. Solids, 51, pp. 22592280, 2003.
 [6] Dao M., Lim C.T. and Suresh S., erratum: Mechanics of the Human Red Blood Cell Deformed with Optical Tweezers, J. Mech. Phys. Solids, 51, pp. 22592280, 2003.
 [7] Eirich F.R., Rheology, theory and applications, vol. 4, Academic Press, 1967.
 [8] Felici M., Filoche M., Sapoval B., Diffusional screening in the human pulmonary acinus, J. Appl. Physiol., vol. 94, pp. 2010, 2003.
 [9] Felici M., Physique du transport diffusif de l’oxygène dans le poumon humain, PhD thesis, Ecole Polytechnique, 2003.
 [10] Fernández M.A., Milisic V., Quarteroni A., Analysis of a geometrical multiscale blood flow model based on the coupling of ODE’s and hyperbolic PDE’s, SIAM J. on Multiscale Model. Simul., vol. 4 (1), pp. 215236, 2005.
 [11] Hsia C.C.W., Chuong C.J.C. and Johnson R.L. Jr, Red cell distortion and conceptual basis of diffusing capacity estimates: finite element analysis, J. Appl. Physiol, vol. 83, pp. 13971404, 1997.
 [12] Hsia C.C.W., Johnson R.L. Jr, Shah D., Red cell distribution and the recruitment of pulmonary diffusing capacity, J. Am. Physio. Soc., vol. 86, pp. 14601467, 1999.
 [13] Jeong J.H. et al, Measurement of RBC deformation and velocity in capillaries in vivo, Microvascular Research, vol. 71, pp 212217, 2006.
 [14] Lenormand G., Henon S. Richert A., Simeon J. and Gallet F., Direct measurement of the area expansion and shear moduli of the human red blood cell membrane skeleton, Biophysical J., vol. 81, pp. 4356, 2001.
 [15] Mauroy B., Hydrodynamique dans les poumons, relations entre flux et gémoétries, PhD thesis, ENS Cachan, 2004. http://www.cmla.enscachan.fr/mauroy/publis/mauroy_these.pdf
 [16] Mills J.P.,Qie L., Dao M., Lim C.T. and Suresh S., Nonlinear Elastic and Viscoelastic Deformation of the Human Red Blood Cell with Optical Tweezers, MCB, vol. 1, no. 3, pp. 169180, 2004.
 [17] Nabors K.L. et al, Red blood cell orientation in pulmonary capillaries and its effect on gas diffusion, J. Appl. Physiol., vol. 94, pp. 16341640, 2003.
 [18] Noguchi H. and Gompper G., Shape transitions of fluid vesicles and red blood cells in capillary flows, PNAS, vol. 102, no. 40, pp. 145914164, 2005.
 [19] Quemada D., Towards a unified model of elastothixotropy of biofluids, Biorheology, vol. 21, pp. 423436, 1984.
 [20] Secomb T.W., Hsu R. and Pries A.R., Motion of red blood cells in a capillary with an endothelial surface layer: effect of flow velocity, AM; J. Physiol. Heart Circ. Physiol, vol. 281, pp. 629636, 2001.
 [21] Tsukada K. et al, Direct Measurement of Erythrocyte Deformability in Diabetes Mellitus with a Transparent Microchannel Capillary Model and HighSpeed Video Camera System, Microvascular Research, vol. 61, pp. 231239, 2001.
 [22] Uzoigwe C., The human erythrocyte has developed the biconcave disc shape to optimise the flow properties of the blood in the large vessesls, Medical Hypotheses, vol. 67, pp. 11591163, 2006.
 [23] Weibel E.R., The Pathway For Oxygen, Harvard University Press, 1984.
 [24] Winslow A.M., Numerical solution of the quasilinear Poisson equation in a nonuniform triangle mesh, Journal of Computational Physics, vol. 135(2), pp. 128138, 1997.
 [25] Yeoh O.H., Characterization of Elastic Properties of CarbonBlackFilled Rubber Vulcanizates, Rubber Chem. Technol., vol. 63, pp. 792805, 1990.