|Home | About | Journals | Submit | Contact Us | Français|
Concentrated erythrocyte (i.e., red blood cell) suspensions flowing in microchannels have been simulated with an immersed-boundary lattice Boltzmann algorithm, to examine the cell layer development process and the effects of cell deformability and aggregation on hemodynamic and hemorheological behaviors. The cells are modeled as two-dimensional deformable biconcave capsules and experimentally measured cell properties have been utilized. The aggregation among cells is modeled by a Morse potential. The flow development process demonstrates how red blood cells migrate away from the boundary toward the channel center, while the suspending plasma fluid is displaced to the cell free layer regions left by the migrating cells. Several important characteristics of microscopic blood flows observed experimentally have been well reproduced in our model, including the cell free layer, blunt velocity profile, changes in apparent viscosity, and the Fahraeus effect. We found that the cell free layer thickness increases with both cell deformability and aggregation strength. Due to the opposing effects of the cell free layer lubrication and the high viscosity of cell-concentrated core, the influence of aggregation is complex but. The lubrication effect appears to dominate, causing the relative apparent viscosity to decrease with aggregation. It appears therefore that the immersed-boundary lattice Boltzmann numerical model may be useful in providing valuable information on microscopic blood flows in various microcirculation situations.
Erythrocytes (i.e., red blood cells, RBCs) are an important determinant of the rheological properties of blood because of their large number density (~ 5 × 106/mm3), particular mechanical properties, and aggregation tendency. Typically, a human RBC has a biconcave shape of ~ 8 μm in diameter and ~ 2 μm in thickness, and is highly deformable [1, 2]. The interior fluid (cytoplasm) has a viscosity of 6 cP, which is about 5 times of that of the suspending plasma (~ 1.2 cP). At low shear stress, RBCs can also aggregate and form one-dimensional stacks-of-coins-like rouleaux or three-dimensional (3D) aggregates [1, 3, 4]. The process is reversible and particularly important in the microcirculation, since such rouleaux or aggregates can dramatically increase effective blood viscosity. RBCs may also exhibit reduced deformability and stronger aggregation in many pathological situations, such as heart disease, hypertension, diabetes, malaria, and sickle cell anemia . The underlying mechanism of RBC aggregation is still under debate [1, 4].
As a consequence of red cell size, the nature of blood flow changes greatly with the vessel diameter. In vessels larger than 200 μm, the blood flow can be accurately modeled as a homogeneous fluid. However, in vessels smaller than this, such as arterioles and venules, and especially capillaries (4–10 μm, i.d.), the RBCs have to be treated as discrete fluid capsules suspended in the plasma. When flowing in microvessels, the flexible RBCs migrate toward the vessel central region, resulting in a RBC concentrated core in the center, and a cell free layer (CFL) is developed near the vessel wall. Such phenomena are important in microcirculation since they are closely related to the flow resistance and biological transport . Traditionally, to incorporate the CFL effects, microscopic blood flows are modeled as multiphase flows with an interface between the peripheral CFL and the RBC concentrated core . Using such an approach, it is difficult to incorporate the effects of cell properties on hemodynamics and it is not suitable for investigation of the microscopic flow structures.
Benefiting from advances in computer and simulation technologies, now it is possible to model individual RBCs as discrete particles suspended in the plasma media. For example, Pozrikidis and coworkers [6–8] have employed the boundary integral method for Stokes flows to investigate RBC deformation and motion in shear and channel flows; Eggleton and Popel  have combined the immersed boundary method (IBM)  with a finite element treatment of the RBC membrane to simulate large 3D RBC deformations in shear flow. Recently, a lattice Boltzmann approach has also been adopted for RBC flows in microvessels, where the RBCs were represented as two-dimensional (2D) rigid particles . Bagchi  has simulated a large RBC population in vessels of size 20–300 μm. However, RBC aggregation, which has been demonstrated experimentally as an important factor in determining hemodynamic and hemorheological behaviors in microcirculation [1, 3, 13], was not considered in these studies. Liu et al.  modeled the intercellular interactions using a Morse potential, thus accounting for RBC aggregation. The cell membrane was modeled with 3D elastic finite elements. Therefore the empirical RBC membrane constitutive relationships, which assume the membrane is a 2D sheet, cannot be utilized directly. Bagchi et al.  have extended the IBM approach of Eggleton and Popel  to two-cell systems and introduced the intercellular interaction using a ligand-receptor binding model. For rigid particles, Chung et al.  have studied the behavior of two rigid elliptical particles in a channel flow utilizing the theoretical formulation of depletion energy proposed by Neu and Meiselman . It has been noted that such a description yields a constant (instead of a decaying) attractive force with large separations, which is not physically realistic . Sun and Munn  have also improved their lattice Boltzmann model by including an interaction force between rigid RBCs.
Recently, we have developed an integrated model to study RBC behaviors in both shear and channel flows [19, 20]. This model combines the lattice Boltzmann method (LBM) for fluid flows and IBM for the fluid-membrane interaction. Crucial factors, including the membrane deformability, viscosity difference across the RBC membrane, and intercellular aggregation, have been carefully considered in the model. The aim of this study was to investigate the effect of different RBC deformabilities and aggregation strengths on CFL structure and flow resistance in microscopic blood flows. We also explored the process of RBC flow development upon the application of an external pressure gradient. Finally, our results are also compared to previous experimental observations and numerical simulations.
The undeformed human RBC has a biconcave shape as described by the following equation :
where c0= 0.207, c1= 2.002, and c2= 1.122. The non-dimensional coordinates (, ) are normalized by the radius of a human RBC (3.91 μm). The RBC membrane is a neo-Hookean, highly deformable viscoelastic material with a finite bending resistance that becomes more profound in regions with large curvature [22, 23]. For the 2D RBC model in the present study, following Bagchi et al. [12, 15], the neo-Hookean elastic component of membrane stress can be expressed as
where Es is the membrane elastic modulus and ε is the stretch ratio. In addition, the bending resistance can also be represented by relating the membrane curvature change to the membrane stress with [7, 15]
where Eb is the bending modulus and κ and κ0 are, respectively, the instantaneous and initial stress-free membrane curvatures. l is a measure of the arc length along the membrane surface. Therefore, the total membrane stress T induced as a result of the cell deformation is a sum of the two terms discussed above:
Here n is the local normal direction on the membrane.
Currently, there are two theoretical descriptions of the aggregation mechanism: the bridging model and the depletion model [1, 4]. The former assumes that macromolecules, such as fibrinogen or dextran, can adhere onto the adjacent RBC surfaces and bridge them together [24–26]. The depletion model attributes the RBC aggregation to a polymer depletion layer between RBC surfaces, which is accompanied by a decrease of the osmotic pressure [27, 28]. Detailed discussions of these two models can be found elsewhere (see, for example, a review by Baumler et al. ). To incorporate the aggregation, here we modeled the intercellular interaction energy as a Morse potential:
where r is the surface separation, r0 and De are, respectively, the zero force separation and surface energy, and β is a scaling factor controlling the interaction decay behavior. The interaction force from such a potential is its negative derivative, i.e., f(r) = −/r. This Morse potential was selected for its simplicity and realistic physical description of the cell-cell interaction; however, it does not imply any particular underlying mechanism. Other aggregation models can be readily incorporated into our algorithm.
A recently proposed IBM-LBM algorithm [19, 20] was employed in this study. In this model, the fluid dynamics was solved by application of LBM while the fluid-membrane interaction was incorporated using IBM. In addition, the fluid property (viscosity here) also has to be updated dynamically with changing RBC shapes and locations. For this purpose, an index θ was introduced, which is 0 in plasma and 1 in cytoplasm. A detailed description of this numerical scheme is provided in the Appendix.
Simulations were carried out over a 100×300 D2Q9 grid (19.5 μm×58.5 μm) domain with 27 biconcave RBCs (Figure 1a). No-slip and periodic boundary conditions were imposed on the channel surfaces and the inlet/outlet, respectively . The volume conservation of RBCs in our 2D simulations was naturally accomplished in IBM. The actual changes in the enclosed RBC areas were typically less than 1% due to numerical errors. Flow was initiated by a pressure gradient along the channel, using a recently proposed boundary treatment . The applied pressure gradient here is 62.5 kPa/m; such a pressure gradient, according to the Poiseuille law, will generate a mean velocity of 1.65 mm/s for a pure plasma flow with no RBCs. This velocity is in the typical range in microcirculation for microvessels of this size [11, 31].
The membrane of each RBC was represented by 100 equal-length segments, and the tube hematocrit (HT) in our simulations was 32.2%.. The simulation parameter values are listed in Table I. Experimental values are not available for the parameters in the Morse potential Eq. (5) for the intercellular interaction. The De values utilized for moderate and strong aggregations were selected based on our previous studies [19, 20]. According to Neu and Meiselman , the zero force distance r0 is ~ 14 μm. In our simulations, r0 was chosen as 2.5 lattice units (0.49 μm) due to the limitation of the grid size and the immersed boundary treatment. To improve the computational efficiency, as is typical in molecular dynamics simulations, a cut-off distance was also adopted. The scaling factor β was adjusted so that the attractive force beyond the cut-off distance rc = 3.52 μm can be neglected.
To visualize the flow field over a time period, we employed virtual tracers in our simulations. The tracer position xt(t) was updated with time t using the Euler scheme
where Δt is the time step, and u is the fluid velocity. Such tracers serve the purpose of flow visualization and have no influence on simulation results. The volumetric flow rates of the blood suspension Qtotal and of the RBCs QRBC through the channel were calculated from the streamwise velocity u and index θ distributions by
respectively. The discharge hematocrit, HD, can also be evaluated from the averaged values of Qtotal and QRBC over a time period when the flow is relatively stable:
In addition, as a measure of the resistance increase due to the RBC existence, the relative apparent viscosity μrel is defined as
where Qp = ΔP H3/12Lμp is the flow rate of a pure plasma flow under the same pressure gradient ΔP/L. Here ΔP is the pressure difference applied across the channel of length L and width H. As we will see in the Results and Discussion section, with the cell migration to the central region, CFLs will be developed near the surfaces. Although this phenomenon have been experimentally observed for many years, the determination of CFL thickness is not trivial, especially in small vessels . Here we define the CFL thickness as the distance from the channel boundary to the θ = 0 position extrapolating downhill from the first two θ ≠ 0 points of the θ profile.
Simulation snapshots at 10 representative time steps are shown in Figure 1 for normal RBCs with moderate aggregation. Once the pressure gradient is applied, the suspension begins to flow, accompanied by large cell deformations due to various local flow conditions. Meanwhile, during the flow development, RBCs near the boundaries migrate toward the channel center. As a result, CFLs are formed near the channel boundaries and the central region has a higher RBC concentration (Figure 1j). Also displayed in Figure 1 is a membrane marker (the black circle) showing the tank-treating-like membrane motion as the cell drifts and deforms.
Six representative fluid tracers, two inside RBCs and four in plasma, are displayed in the snapshots in Figure 1. The membrane marker (black circle) and the fluid tracer inside that RBC undergo a tank-treading motion. The four tracers in plasma are labeled as 1 to 4. For tracer 1 initially positioned near the lower boundary, its movement is mainly in the flow direction with small transverse fluctuations. However, tracers 2 and 3, originally between RBCs, are quickly displaced toward the boundary CFL regions by the center-ward migrating RBCs (see Figures 1b–g for tracer 2 and 1e–i for tracer 3). Once the tracers arrive in the CFL regions, they behave very similarly to tracer 1 and move mainly in the flow direction. Tracer 4 is trapped between RBCs in the center region and has never entered the CFL regions during the simulation. Their trajectories are also plotted in Figure 2 to show their different behaviors more clearly. As expected, tracer 4 moves much faster than tracers 1–3 at the periphery., In general, a fluid particle will be displaced to the CFL region by migrating RBCs, unless it is closely surrounded by RBCs. Once a fluid particle reaches the CFL region, it will generally maintain that transverse location as it flows downstream. Although this analysis involves virtual fluid tracers, such results imply that small solid particles in blood will accumulate in CFL regions, as has been observed in experiments with platelets .
A more quantitative analysis for the developing process discussed above can also be performed. Figure 3a displays the volumetric flow rates of the blood suspension Qtotal and of the RBCs QRBC through the channel. When the pressure gradient is applied, the blood suspension begins to flow forward and both flow rates increase gradually. Unlike that for a pure fluid, the total flow rate Qtotal exhibits random fluctuations and only a quasi-steady stable state can be established, since the RBC configuration is constantly changing. Such fluctuations are more profound in the RBC flow rate QRBC , because any RBC configuration change will affect both the velocity and index distributions, see Eq. (8). The resulting discharge hematocrit value for this simulated case, HD = 38.9%, is higher than the tube hematocrit HT = 32.2%, and consistent with the well-known Fahraeus effect.
Figure 3b shows the RBC distribution profiles during the flow development, where θ has been averaged along the channel. It can be seen that, as time progresses, a layer of θ = 0, which by definition is pure plasma, develops near each boundary surface and its thickness increases with time. Also displayed in Figure 3c are the streamwise averaged velocity profiles at several time points. As observed in in-vivo and in-vitro experiments , such velocity profiles have blunt shapes due to the high fluid viscosity of the RBC-rich core in the central region. Compared to a pure plasma flow under the same pressure gradient (Figure 3c, dashed line), the flow rate of blood suspension Qtotal is smaller, and the relative apparent viscosity μrel = 1.29 is larger than unit, indicating an increase in flow resistance due to the RBC presence.
Similar 2D systems have been investigated by Sun and Munn  and Bagchi . Aggregation was not considered in those studies. For rigid RBCs with HT = 30.5% flowing in a 20 μm channel, Sun and Munn  found that the relative apparent viscosity μrel was 1.50 and hematocrit ratio HT /HD was 0.90. Such values are close to our results for less deformable RBCs with no aggregation, i.e., μrel = 1.66 and HT /HD = 0.875. The differences may be due to the different tube hematocrit and RBC shapes, as well as the slight deformability in our simulation. Also no clear CFLs were observed in Ref.  as in the top-right panel in our Figure 4. Recently, Bagchi  employed a similar RBC model and studied deformable RBCs with HT = 33% also in a 20 μm channel. The CFL thickness was about 2.8 μm, which is in a good agreement with the two-phase model by Sharan and Popel , and is close to our estimate of δ = 2.44 μm for the normal RBCs with no aggregation. However, the relative apparent viscosity μrel = 1.60 obtained by Bagchi is much larger than our result μrel = 1.27 The reason for this discrepancy is not clear One possibility is the different mean flow viscosities as discussed by Bagchi.
Unlike these previous 2D numerical studies, we have not attempted to compare our simulation results to the empirical relationships proposed by Pries et al. . We believe that 2D studies, although valuable, can only provide qualitative information. For example, the CFL/RBC-core cross-sectional area ratio are 2δ/H and 4δ (D − δ)/D2 (D is the vessel diameter) for 2D and 3D situations, respectively. For δ = 3 μm and H = D = 20 μm, 2δ/H = 0.3 and 4δ (D − δ)/D2 = 0.51, implying a less profound CFL lubrication effect and a lower core hematocrit in 2D simulations compared to the corresponding 3D systems. Also the apparent viscosity of a RBC suspension may be very different in 2D and 3D systems, even with the same hematocrit. It is interesting to note that Sun and Munn  found a good agreement to the Pries formula  for the hematocrit-viscosity relationship, although the cells were represented as rigid particles. However, for HT = 30.5% in a 20 μm channel, the simulated discharge hematocrit HD by Sun and Munn was 33.9%, while the prediction according to Pries et al.  is 42.5%. In Ref.  Bagchi reported a similar agreement of the HD ~ μrel relationship to the Pries empirical formula. However, the discharge hematocrit HD was not obtained from simulations, although it is readily available from the velocity and index distributions as shown in Eq. (9).
The simulations with different RBC deformabilities and aggregation strengths (Table I) are summarized in Figures 4 and and55 and Table II. Less deformable RBCs are modeled by increasing the membrane shear and bending modular to twenty times of their normal values, i.e., Es = 1.2 × 10−4 N/m and Eb = 4.0 × 10−18 Nm. Three different aggregation strengths, De = 0 (none), 5.2 × 10−8 (moderate), and 1.3 × 10−6 μJ/μm2 (strong), were employed based on previous numerical studies . In Figure 4, we plot the quasi-steady states of the six situations with different deformabilities and aggregation strengths. Clearly, many less deformable RBCs (right column in Figure 4) preserve their original biconcave shape and the deformations of others are less evident, when compared with their counterparts with normal RBCs (left column in Figure 4). Another feature is the presence of CFLs near the channel boundaries: Except for the case with less deformable RBCs and no aggregation capability (Figure 4, top-right), obvious CFLs can be seen. The CFL thickness increases with both the aggregation and the cell deformability, since a stronger aggregability tends to attract RBCs closer and generate a more concentrated RBC core. For less deformable RBCs, the membrane rigidity prevents large deformation and large cell-cell contact area, resulting in a less dense RBC core and thinner CFLs than with normal RBCs. Similar behaviors have been reported in previous studies [11, 12]. Also displayed in Figure 4 are the blunt velocity distributions.
In Figure 5, we have plotted the averaged fluid index θ and streamwise velocity u for these six different cases. CFLs (θ = 0 regions in Figure 5a) have developed near the boundaries and concentrated RBC cores with high θ values have formed in the central region. The measured CFL thicknesses δ are listed in Table II. Clearly δ increases with both the cell deformability and aggregation strength. For the less deformable RBCs with no aggregation, no CFLs are evident in Figure 4b, similar to the findings of Sun and Munn , where RBCs were modeled as rigid particles. The measured δ for this case is 1.32 μm, the smallest value in all six cases. The relative apparent viscosities calculated according to Eq. (10) are listed in Table II. All these relative viscosities are greater than 1 because of the existence of RBCs. Obviously the RBC deformability plays a definitive role. The relative apparent viscosity increases with a decrease in cell deformability. However, the aggregation effect is more complex. With stronger aggregation, the denser RBC core has a higher viscosity; whereas the larger CFL will generate a more significant lubricating effect, which tends to reduce the global apparent viscosity. This finding may help to explain the contradictory experimental findings in the literature . In our simulations at low flow rates, the lubrication effect appears more profound, and the relative apparent viscosity μrel generally decreases with aggregation. This trend is in a good agreement with experimental findings by several groups [36–39]. However, when moderate aggregation is introduced, μrel indeed increases slightly compared to the case with no aggregation, indicating that the larger core viscosity has suppressed the lubrication benefit from thicker CLFs. We have also calculated the discharge hematocrit HD by Eq. (9) and the tube-discharge hematocrit ratio HT /HD. The tube hematocrit HT is the same, 32.3%, for all cases, and the discharge hematocrit HD varies from 36.8% to 40.0%. The corresponding HT /HD ratios changes from 0.875 to 0.805, all less than 1, the well-known Fahraeus effect.
A recently developed IBM-LBM model has been applied to study the effects of RBC deformability and aggregation on the hemodynamic and hemorheological behaviors of microscopic blood flows. When a pressure gradient is applied, RBCs near the channel boundaries migrate toward the channel center, creating a concentrated RBC core, and plasma fluid in that region is displaced toward the CFL regions. Tank-treading motion of the RBCs also occurred. We found that a fluid particle is likely to be displaced laterally by RBCs, unless it is closely surrounded by RBCs. This finding may help to explain the typical U-shaped distribution of platelets in microvessels observed experimentally . In addition, several typical microscopic hemodynamic and hemorheological characteristics have been observed, such as the blunt velocity profile, changes in apparent viscosity, and the Fahraeus effect. We also found that CFL thickness increases with cell deformability, since the membrane deformability allows the cells to deform and create a concentrated core. However, the influence of aggregation is more complex due to the opposing effects of the CFL lubrication and the high RBC core viscosity. This leads to a decrease in effective viscosity with higher levels of aggregability in our study, in agreement with in vitro experimental findings . Overall, we find that our model predicts important aspects of blood flow behaviors observed both in microcirculatory vessels and in vitro. This suggests that the model may also be helpful in understanding the factors involved, for example, in the intravascular distribution and delivery to target organs of stem cells, drugs and genetic materials that are less amenable to experimental study. Extending this model to 3D systems would provide more quantitative results.
This work was supported by NIH grants HL18292 and HL079087. J. Zhang acknowledges the financial support from the Natural Science and Engineering Research Council of Canada (NSERC) and a start-up grant from the Laurentian University. P. C. Johnson is also a Senior Scientist at La Jolla Bioengineering Institute. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET:www.sharcnet.ca).
In LBM, a fluid is modeled as pseudo-particles moving over a lattice domain at discrete time steps. The major variable in LBM is the density distribution fi(x, t), indicating the particle amount moving in the i-th lattice direction at position x and time t. The time evolution of density distributions is governed by the so-called lattice Boltzmann equation, which is a discrete version of the Boltzmann equation in classical statistical physics [29, 40]:
where ci denotes the i-th lattice velocity, Δt is the time step, and Ωi is the collision operator incorporating the change in fi due to the particle collisions. The collision operator is typically simplified by the BGK  single-time-relaxation approximation:
where τ is a relaxation parameter and the equilibrium distribution can be expressed as i 
Here ρ = Σifi is the fluid density and u = Σifi ci/ρ is the fluid velocity. Other parameters, including the lattice sound speed cs and ti, are lattice structure dependent. Through the Chapman-Enskog expansion, one can recover the macroscopic continuity and momentum (Navier-Stokes) equations from the above-defined LBM microdynamics:
where ν is the kinematic shear viscosity given by
and P is the pressure expressed as
It can be seen from the above description that the microscopic LBM dynamics is local (i.e., only the very neighboring lattice nodes are involved to update the density distribution fi), and hence a LBM algorithm is advantageous for parallel computations . Also its particulate nature allows it to be relatively easily applied to systems with complex boundaries, such as flows in porous materials . Such features could be valuable for future studies of RBC flows in more complicated and realistic systems.
In the 1970s, Peskin  developed the novel immersed boundary method (IBM) to simulate flexible membranes in fluid flows. The membrane-fluid interaction is accomplished by distributing membrane forces as local fluid forces and updating membrane configuration according to local flow velocity.
Figure A.1 displays a segment of membrane and the nearby fluid domain, where filled circles represent membrane nodes and open circles represent fluid nodes. In IBM, the membrane force f(xm) at a membrane marker xm induced by membrane deformation is distributed to the nearby fluid grid points xf through a local fluid force
Here h is the lattice unit. The membrane velocity u(xm) can be updated in a similar way according to the local flow field:
We note that both the force distribution Eq.(A.8) and velocity interpolation Eq.(A.11) should be carried out in a 4h × 4h region (dashed square in Figure A.1) , instead of a circular region of radius 2h as described elsewhere [43–45]. This is due to the specific approximation of the delta function Eq.(A.9) adopted in IBM. It can be shown that the sum of non-zero values of Eq.(A.9) in the square region is 1, and hence missing any nodes inside the square and outside of the circle of radius 2h would produce an inaccuracy in fluid forces and membrane velocity.
The fluids separated by the membrane can have different properties, such as density and viscosity. In our system, the densities of cytoplasm of RBC and plasma are similar and the difference could be neglected for low Reynolds numbers. However, for normal RBCs, the viscosity of cytoplasm is 5 times of that of the suspending plasma. As the membrane moves with the fluid flow, the fluid properties (here the viscosity) should be adjusted accordingly. This process seems straightforward; however, it is not trivial. Typically, an index field is introduced to identify the relative position of a fluid node to the membrane. In previous IBM studies, Tryggvason et al.  solved the index field through a Poisson equation over the entire domain at each time step, but the available information from the explicitly tracked membrane is not fully employed . On the other hand, Shyy and coworkers [43–45] suggested to update the fluid properties directly according to the normal distance to the membrane surface. We prefer the latter approach for its simple algorithm and computational efficiency. However, the originally proposed formulations of the Heaviside function are coordinate dependent . To overcome this drawback, here we modify this approach by simply defining an index d of a fluid node as the shortest distance to the membrane. If the node is close to two or more membrane segments, the index taken is that to the closest one. This index also has a sign, which is negative outside of the membrane and positive inside. Similar to the approach of Shyy and coworkers [43, 44], a Heaviside function is defined as
A varying fluid property α can then be related to the index d by
where the subscripts in and out indicate, respectively, the bulk values inside and outside of the membrane. It can be seen from Eq.(A.12) that, as the membrane moves, only the nodes close to the membrane (|d| ≤ 2h) will be affected. Therefore, in the computer program we need only to calculate the index and update the fluid properties in the vicinity of the membrane. Such a process should be more efficient than solving a Poisson equation over the entire fluid domain .
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.