The process of development of an embryo and subsequent functioning of the cardiovascular system is important for the explanation of many phenomena occuring within this system. Vertebrates have a functional vasculature, with a heart that pumps blood and blood vessels that have a clearly defined endothelium. However, observation of the blood vessels in living embryos is difficult either because the embryos are developing within the uterus of the mother or because they are not transparent enough. A fish species called teleostei, more precisely a subspecies called zebrafish (in Latin Danio rerio), has great advantages over other species for studying vascular development. First of all, the dimension of zebrafish is very small, since an adult zebrafish can grow to 4 cm on average. Zebrafish are very fertile and lay a large number of eggs. The reproduction is external, outside the mother's body. Optical analysis of this species is also very easy. Due to its small dimensions, the embryos can survive with only a small amount of oxygen that they receive by passive diffusion. All the interior organs, such as eyes, brain, heart, inner ear are developed within the first 3 days post fertilization. The cardiovascular system is one of the first systems that is formed during the embryonic development of zebrafish. After only 1 day post fertilization, the zebrafish is sufficiently developed for the blood to start circulating. Detailed explanation of the formation of all major blood vessels in zebrafish can be found in the literature (Isogai et al., 2001; Ellertsdóttir et al., 2010). Even though there are certain variations in details of the anatomy of the cardiovascular system of this fish, the basics of the system are similar to other vertebrates. Because of all the mentioned advantages, many papers in the literature analyze many aspects of development of the cardiovascular system of zebrafish, and present numerous genetic analysis of mutation of diverse genes (Stainier et al., 1995; Weinstein et al., 1995). The zebrafish species has also been used for many diverse investigations, e.g., as a model to study angiogenesis (Chávez et al., 2016) or to examine the pathophysiology of myofibrillogenesis and muscular dystrophies (Raeker et al., 2014) etc. The cardiac and metabolic physiology of zebrafish was analyzed in the literature (Gore and Burggren, 2012), as well as development of the zebrafish heart in vivo (Hou et al., 2014).

There are also numerous studies that study the formation of erythrocytes in vertebrates, as well as in teleostei (Swaen and Brachet, 1899; Strawinski, 1949; Vernier, 1969). The process of generating blood cells is known as hematopoiesis. Zebrafish erythropoiesis begins in the mesodermal layer during embryonic development (Kulkeaw and Sugiyama, 2012). Zebrafish hematopoiesis undergoes two waves: primitive and definitive. Between 12 and 24 h post fertilization, primitive hematopoiesis starts in the intermediate cell mass (ICM), located between the somites and yolk sac. During this primitive wave, erythrocytes and macrophages are produced. By 24 h post fertilization, primitive erythrocytes enter circulation and mature. After maturation the erythrocytes retain the nucleus, elliptical shape and express hemoglobin (de Jong and Zon, 2005; Li et al., 2014).

The definitive hematopoiesis generates hematopoietic stem cells (HSCs) that differentiate into erythrocytes, lymphocytes, and platelets. It is also called adult hematopoiesis and it has the capability of self-renewal and produces all mature hematopoietic lineages (Falenta and Rodaway, 2011; Jin and Wen, 2011). In zebrafish, the Runx1 transcription factor is crucial for HSCs formation, which is observed at early 24 h post fertilization. The Gata-1, transcription factor expression is vital for primitive hematopoiesis and it is found mainly in the lateral plate mesoderm that migrates medially during the formation of ICM. Gata-1 positive cells are expressed during the differentiation of ICM to proerythroblasts (Li et al., 2014).

Using confocal microangiography (Weinstein et al., 1995) or green fluorescent protein (GFP, Motoike et al., 2000; Lawson and Weinstein, 2002), and relying on the fact that zebrafish are small and very transparent, it is possible to image the blood vessels throughout the entire depth of the zebrafish. This way the anatomy of the blood vessels is obtained and the flow of red blood cells (RBCs) through the blood vessel is monitored. Unlike some other techniques used for the analysis of the vasculature, the mentioned approach does not jeopardize in any way the fish and the obtained images are related to the fully active blood circulation.

Due to the important role of RBCs in vertebrates and human organisms, the analysis of the dynamics of motion of these cells separately is one of the most important problems in physiology and biomechanics. Many authors have analyzed the behavior of RBCs, both experimentally and theoretically. Behavior of synthetic capsules has been experimentally observed (Chang and Olbricht, 1993; Walter et al., 2000) and similar experiments were performed with RBCs (Gaehtgens et al., 1980; Pries and Secomb, 2011). Many authors have investigated the mechanical properties of erythrocytes, with a special focus on the characteristics of the cellular membrane, in the past century (Skalak, 1976; Hochmuth and Waugh, 1987), as well as in the past decade (Mukhopadhyay et al., 2002; Kuzman et al., 2004; Li et al., 2005). A theoretical model was used to simulate the motion of RBCs through capillaries with variable cross-sections, in order to predict the resistance of the vessel to the motion and deformation of RBCs in living microvessels (Secomb and Hsu, 1996). Recently, numerical simulations were performed on idealized arteriole-sized blood vessels and the influence of motion of RBCs on shear stress on the blood vessel walls was analyzed (Gambaruto, 2016). Secomb et al. (2007) performed experiments to analyze the behavior of human RBCs within a glass tube, that has a diameter smaller than the diameter of the RBC. This is similar to the conditions in human capillary blood vessels, which are also narrower than a single RBC. In this study, they observed the capability of an RBC to adapt to the changes in geometry during the flow and to change its shape significantly. Numerical simulations were also performed (Djukic and Filipovic, 2015) and compared with experimental results presented by Secomb et al. (2007). Through this comparison it was demonstrated that the proposed numerical model is capable of accurately predicting the change of shape that the RBC undergoes during its motion through a narrow glass tube. Most of the numerical models previously published in the literature consider only the motion of particles and RBCs in simpler geometrical conditions, where the high defomability of RBCs does not come to the fore. The goal of this paper is to model the motion of RBCs through a complex geometrical domain.

In this paper the behavior of RBCs within the caudal vein plexus during embryonic development of zebrafish is analyzed. The embryonic capillary plexus has its honeycomb-like appearance due to the aggregation of many transluminal intussusceptive pillars. The blood flow videos of observed phenomena in the capillary plexus regions were tracked in vivo, to isolate individual RBC movement. Subsequently, numerical simulations were performed in geometries and under conditions that are defined according to the experimental setup. The results are used to predict and analyze the motion and deformation of RBCs.

The paper is organized as follows: Section Materials and Methods describes the methods that were used to obtain experimental data and the numerical model that was used in numerical simulations. Section Results presents the results obtained in numerical simulations and the comparison with experimental data. Numerical methods are discussed in Section Discussion and Conclusion, including the conclusions about the significance of the presented results and applications of the proposed numerical model.

In this paper it is considered that the RBC of the zebrafish has an elliptic shape, according to the previous investigations of the zebrafish cardiovascular system (Watkins et al., 2012). The model of the RBC membrane that is discretized into a triangular mesh and that is used in simulations presented in this paper is shown in Figure 2. The cross-section of RBC along the plane that is parallel to x0z plane and that contains the center of gravity of the RBC is also shown in Figure 2.

Because of the structure of RBCs and all the conditions to which they are exposed during its motion through the circulation, four types of resistance to deformation are defined. This was discussed in Section Modeling the Deformation of the RBC. Values of the parameters that define these reactions to deformation have been measured in the literature (Dao et al., 2003; Bagchi et al., 2005). For modeling the behavior of RBCs in this paper, the values are taken to be equal to the ones proposed by Dupin et al. (2007).

The parameter of resistance to the change of volume is equal to K^V = 50pNμm⁻¹. The parameter of resistance to the local change of membrane area is equal to K^Al = 1.67·10⁻¹pNμm⁻¹, while the parameter of resistance to the global change of membrane area is equal to K^At = 1.67pNμm⁻¹. The parameters that were used to define the Skalak material model of cellular membrane are equal to ks=0.75·103pNμm-1 and kα=75·103pNμm-1, according to the literature (Skalak et al., 1973). The parameter of resistance to the bending of the membrane is equal to K^B = 10⁻¹pNμm⁻¹.

Simulations using the proposed numerical model were performed for two cases of motion of an individual RBC in the caudal vein plexus of a living zebrafish whose circulation was observed 32 h post fertilization, for an overall time period of 11 s. The geometry of the fluid domain is created directly from the microscopic experimental images. The width of the vein plexus is equal to 100μm. The bounding walls of the fluid domain (the vessel walls) and the intussusceptive pillars are modeled using the Bounce-back approach (Ginzbourg and d'Humières, 1996; Gallivan et al., 1997). This practically means that all particles that collide with the walls, return to the fluid domain with the same velocity, which imposes that the velocity at the walls is equal to zero. In the observed part of the capillary plexus blood is flowing from left to right. The outflow boundary condition is defined at the right boundary wall. This means that normal derivatives at the boundary of the relevant quantities are set to be equal to zero. The left boundary wall of the simulation domain is used to define the inlet velocity of the blood. The velocity profile at the inlet at the beginning of the simulation is defined according to the profile in Poiseuille flow. The value of velocity that is prescribed at nodes that are located on the left boundary of the lattice mesh is calculated using experimental data. The blood flow videos of the capillary plexus region were used to track motion of individual RBCs using imaging techniques. By analyzing the motion of these RBCs, their velocity is calculated and this value is used to define the velocity of the blood flow at the inlet. The maximum value of velocity at the inlet that was used in numerical simulations was equal to 60μm/s.

Figure 3 shows the results for the first considered initial position of the RBC, for several moments in time. On the left the microscopic images of the caudal vein plexus are shown, whereas the considered RBC is denoted by a red line, highlighted in a blue circle. Yellow dotted lines mark the boundaries of the vessel and the pillars within the vessel. On the right the results obtained using numerical simulations are shown. In the middle only the shapes of the considered RBC are isolated for easier comparison (the shape obtained in numerical simulation is denoted by blue color and the shape obtained from the microscopic image is denoted by red color). The case of motion of the RBC that is shown in Figure 3, exhibits the attachment of the RBC to the vascular wall and its subsequent rolling. Figure 4 shows the comparison of results, for the second considered initial position of the RBC, for several moments in time. In this case, the bending and hanging of the RBC on the vascular pillars as well as flowing away and restoration of the original form of the RBC are shown. The shapes obtained in experiments and numerical simulations are also quantitatively compared. The values of area of the cross-section of the RBC along the plane that is parallel to the direction of motion and that contains the center of gravity of the RBC are compared. Table 1 shows the percentage error, that was obtained in numerical simulations, compared to the values that were calculated using experimental data, for the first initial position of the RBC. Table 2 shows the percentage error of the calculated areas of the cross-section of the RBC that was obtained in numerical simulations, compared to the values that were calculated using experimental data, for the second initial position of the RBC. As it can be seen from the isolated shapes of the considered RBCs and the values of errors given in Tables 1, 2, the results obtained using numerical simulation agree well with experimental results.

Several approaches that deal with the numerical modeling of the behavior of RBCs in the fluid domain are presented in the literature. Pozrikidis (1995) applied the Boundary Element Method (BEM) to model the motion of RBCs. This method was also used by other authors (Kraus et al., 1996; Lac et al., 2004). Ramanujan and Pozrikidis (1998) extended this approach and considered various material models to define the relationship between the reaction force and deformation of the particle. Eggleton and Popel (1998) used the Immersed Boundary Method (IBM) to model the influence of solid on fluid and vice versa. Sui et al. (2008) improved this approach by adding the refinement of the mesh, in order to be able to model the motion more precisely, with higher mesh density near the particle. Discrete particle methods have also been applied to model this type of phenomena. Boryczko et al. (2003) modeled the solid using classic continuum mechanics, while the fluid was modeled as a cluster of particles. On the other hand, based on research published in the literature (Koshizuka et al., 1995), Tsubota et al. (2006) modeled both fluid and solid as a cluster of particles, analyzed their mutual interactions and subsequently modeled the motion of RBCs through the blood plasma. Simulations of a larger number of RBCs and their mutual interactions have been analyzed in several papers (Liu et al., 2004; Dupin et al., 2007; Doddi and Bagchi, 2009). However, in most of the mentioned papers, only the motion of particles and erythrocytes in simpler geometrical conditions has been considered. Under these conditions the feature of erythrocytes to drastically change its shape does not come to the fore. The goal of this paper was to present a model that simulates the motion of RBCs through complex geometrical domains, more precisely through real physical domains obtained experimentally by recording the blood flow through the caudal vein plexus of living zebrafish.

Two initial positions of the RBC within the caudal vein plexus were analyzed. The shapes of the RBC during its motion obtained in numerical simulations agree well with shapes that were extracted from experimental data. The standard deviation of numerical values obtained from the experimental values of the area of the observed cross-section of the RBC is equal to 6.91% for the first considered initial position of the RBC. For the second considered initial position, the standard deviation is slightly greater and is equal to 11.1%. As it can be observed from Tables 1, 2, most values of the obtained error were less than 5%. The greatest error was obtained for the last two moments in time, for the second initial position of the RBC (last two rows in Table 2). These higher values are caused by the fact that in the numerical model a relatively simple approach was used to simulate the intussusceptive pillars. Namely, the pillars were treated as boundary walls, and the RBC was bounced-back from the pillars, just like the fluid particles. If a more detailed interaction between the pillar walls and the RBC was implemented, then this error would be additionally reduced. This will be the main direction of future improvements of the numerical model.

The modeling of oxygen transport and metabolism in organs is significant for the understanding of the functioning of organs and cellular functions from the physiological aspect and the understanding of ischemic and hypoxic conditions. The process of oxygen transport consists of several stages: passage of oxygen molecules through the membrane of the RBC, motion of oxygen together with the RBC through the blood plasma, passage through the vascular wall and arrival to the mitochondria. As stated in the literature, the delivery of oxygen to the tissue is determined mainly by three factors: capillary blood flow, hematocrit, and arterial pressure of oxygen (Li et al., 1997; Dash and Bassingthwaighte, 2006). However, oxygen usage on a cellular or tissue level cannot be measured directly. In order to quantitatively explain phenomena that occur during oxygen transport, it is necessary to use mathematical models and numerical simulations, which have been quite successful in the study of many complex biological systems (Wolkenhauer, 2013).

There have been several methods proposed in literature that mathematically describe the transport of oxygen. Dash and Bassingthwaighte (2006) concluded in their study that a significant decrease in blood flow can cause acute hypoxia and prevent cells from functioning normally. Beyer et al. (2002) proposed a convection-diffusion-reaction mathematical model that simulates the transport of oxygen from RBCs to mitochondria, where they treated RBCs and blood plasma as two separate flows. Li et al. (1997) analyzed oxygen transport at the regional level with imaging techniques using tracer ¹⁵O-oxygen for positron emission tomography. They observed that with a decrease of capillary blood flow, oxygen delivery by flow to the tissue can become inadequate. Also, when they analyzed tracer kinetics, the initial statement was confirmed—when the capillary blood flow or hematocrit or arterial pressure of oxygen were reduced, the retention of tracer oxygen was prolonged, and the fraction of tracer oxygen in the outflow of the observed blood vessel was reduced.

A more extensive model of O₂/CO₂ transport and exchange in the microcirculation, that covers all aspects and all phenomena that happen during this process has yet to be developed. In all previously mentioned models, blood is assumed to consist of two continuous coexisting phases and RBCs are only one “layer” within the blood. The entire microvessel is most commonly described as a tube, where the RBCs are assumed to be located within the central column, while the surrounding region represents blood plasma. The velocity of RBCs and blood plasma can be different within the observed domain and this way the motion of RBCs relative to the blood plasma is taken into consideration. On the other hand, because oxygen is mostly transported through the body within the RBCs, in the form of oxyhemoglobin, the time taken for the transmitted nonextracted oxygen to arrive to the targeted location mainly depends on the velocity of the RBCs. Thus, the aspect of motion of individual RBCs within the blood vessel has to be considered. The method proposed in this paper can be connected to the previous methods modeling the transport of oxygen and can give valuable insights about the motion of RBCs through complex domains of living microvessels.

The numerical model presented in this paper can also be used to analyze erythroid diseases, where disruption of RBC morphology and mechanics occurs. These diseases include hereditary spherocytosis, hereditary elliptocytosis, sickle cell disease etc. The changes in RBCs in these cases have been studied in the literature (Diez-Silva et al., 2010; Fisseha and Katiyar, 2012; Du et al., 2015) and the presented numerical model can be used in combination with these findings to provide additional insight into the mentioned phenomena. The numerical model presented in this paper is capable of simulating the microcirculation, including the complete simulation of dynamics of motion of individual cells. Most numerical models that have been presented in literature so far are limited to simulations in geometrically simple blood vessels. Due to the accuracy of the obtained solutions that was demonstrated through the comparison of the numerical results with experimental results, the proposed model and the developed software can be used for simulations of the complex motion of viscoelastic bodies, such as RBCs, in capillaries with complex geometry.

TD conceived the study, selected data for numerical simulations and developed the numerical model; SK perfomed experiments with zebrafish; TD and IS performed numerical simulations; SK and IS participated in study conception and assisted in writing the manuscript; VD and NF participated in the study coordination, gave their important intellectual content, participated in data interpretation and contributed to the revision of the manuscript. All authors have read and approved the final manuscript.

This paper is supported by grants from Ministry of Education, Science, and Technological Development of the Republic of Serbia (projects number III41007 and ON174028). This paper is also supported by the Swiss National Science Foundation, SCOPES project (JRP: IZ73Z0_152454/1).