Extended photoacoustic transport model for characterization of red blood cell morphology in microchannel flow

: The dynamic response behavior of red blood cells holds the key to understanding red blood cell related diseases. In this regard, an understanding of the physiological functions of erythrocytes is signiﬁcant before focusing on red blood cell aggregation in the microcirculatory system. In this work, we present a theoretical model for a photoacoustic signal that occurs when deformed red blood cells pass through a microﬂuidic channel. Using a Green’s function approach, the photoacoustic pressure wave is obtained analytically by solving a combined Navier-Stokes and photoacoustic equation system. The photoacoustic wave expression includes determinant parameters for the cell deformability such as plasma viscosity, density, and red blood cell aggregation, as well as involving laser parameters such as beamwidth, pulse duration, and repetition rate. The eﬀects of aggregation on blood rheology are also investigated. The results presented by this study show good agreements with the experimental ones in the literature. The comprehensive analytical solution of the extended photoacoustic transport model including a modiﬁed Morse type potential function sheds light on the dynamics of aggregate formation and demonstrates that the proﬁle of a photoacoustic pressure wave has the potential for detecting and characterizing red blood cell aggregation. inﬂuence relatively small in For an application, we studied the red blood cell aggregation in a microcirculatory system. Our detailed analytical solution of the extended photoacoustic transport model combined with a modiﬁed Morse type potential function accounts for the dynamics of aggregate formation. Hence, the distinctive features of the photoacoustic waves enable to understand diseased or damaged RBCs having abnormal shapes and sizes, which can result from a variety of diseases. Our detailed solution of PA wave to detect RBC aggregation is obtained by solving the combined Navier-Stokes and bio-heat transfer equations. This study also presents the spectral features of the PA signals by considering the interaction between RBCs modeled by the Morse potential. We think that our work can be helpful to improve hybrid simulation approach for not only the characterization of red blood cell morphology but also for cancer cell lines in the microchannel ﬂow. To describe and analyze the mechanical processes inside a microﬂuidic device, our hybrid model can pave the way toward the development of new physical models in the rapidly advancing ﬁeld of lab-on-chip technology and photoacoustics. As a future work, in order to analyze the inﬂuence of the deformation of the membrane on pressure wave, we will implement the spring network approach beside the Morse potential.


Introduction
Red blood cells (RBCs) have the deformation capability in microcapillaries. RBC membrane allows cells to squeeze, deform, and reform through capillaries. The geometric shape of RBC provides an efficient transport through narrow capillaries and increases the surface area in order to maximize oxygen transport. The alternation of the cell membrane deformability under some disease conditions is significantly different than healthy ones [1,2]. To illustrate, the RBCs of sickle cell disease have the most significant loss of deformability because of having different morphologies depending on its density [3][4][5][6]. The small changes in the morphology of RBCs are characterized by the absorption-based contrast of photoacoustic (PA) imaging which is complementary to other imaging modalities in terms of contrast mechanism, penetration, spatial resolution, and temporal resolution [7][8][9][10][11][12][13][14][15][16][17][18].
Recent studies seek to provide insights to measure RBC deformability via experimental techniques such as atomic force microscopy (AFM) [19][20][21][22], micropipette aspiration technique [23][24][25][26], optical and magnetic tweezers [27][28][29][30][31][32], microfluidic device techniques [2,[33][34][35]. Moreover, mathematical (analytical and numerical) models have been also presented in order to distinguish cell types and identify several diseases [28,32,[36][37][38]. Using these experimental techniques and mathematical methods, intrinsic mechanical properties of RBCs (geometry and cytoplasmic viscosity), and viscoelastic properties of their membrane cortex structure are assessed. The shear modulus of RBC is inferred using optical tweezers while two microspheres attached to the opposite sides of RBC are trapped. Unless the laser beam is tightly focused on microbeads attached to the cells, the cell damages occur because of temperature rise [39]. The stiffness of cells is measured via AFM by isolating the AFM system from surrounding vibrations [40]. Micropipette aspiration is another technique for investigating cellular mechanical properties via aspirating pressure. Although this technique provides a broad range of applications, the unstable base pressure coming from vapor evaporation in the aspiration chamber may occur owing to the dependence of operator's skills [26]. Furthermore, magnetic tweezer has been applied to monitor mechanical properties of single molecules, intermolecular bonds, and whole cells. It gives information about complex modulus of elasticity and the local viscoelasticity of the cells. However, it is required to ensure highly localized magnetic fields, field gradients, and spherical magnetic beads with various sizes to obtain more accurate viscoelastic properties of the cells [41].
The hemoglobin within RBCs as an endogenous contrast agent absorbs visible and near-infrared light, which enables to image RBCs without staining via PA microscopy [7,9,42,43]. Even though the image does not give information about their shapes and structures, PA signal with amplitude, width, rise, and relaxation time enables to observe morphological changes in RBCs. The changes in the PA response depending on the cell elastic properties differentiate between healthy RBCs and diseased or damaged RBCs [15,44,45]. In the context of the detection of individual RBCs, there are several recent PA studies utilizing theoretical, experimental, and computational methods [46]. Saha et al. firstly described the PA pressure field generated by a collection of erythrocytes with the help of a theoretical model developed by Diebold et al. [46][47][48][49]. In vivo RBC age or surface charge density was presented using a two dimensional simulation study based on the PA signal properties of erythrocyte aggregation level which is related to the suspending medium [47,49]. In another work, Saha and Kolios not only presented experimental results but also showed that the theoretical simulation study assessing PA signals could be followed to determine the degree of aggregation [50]. This frequency domain theoretical model was also used for a computer simulation study in order to differentiate intraerythrocytic stages of malarial parasite such as normal, ring, trophozoite, and schizont stages [46].
Mathematical models to calculate PA signal produced from one or more cells have been reported in a limited number of studies [48, 50-52]. Generally, cells were approximated as fluid spheres. The PA pressure wave produced by a collection of absorbing spheres was usually considered as a linear superposition of spherical waves emitted by individual sources [50]. It was assumed that acoustic waves coming from an absorbing sphere were not influenced by other particles in the medium, and all double and multiple scatterings of light beam were also ignored. The other study utilizing mathematical models was the description of PA signal in the time domain produced from one or more cells using an absorption model. The results of this study demonstrated that PA amplitude does not depend on the position of the cells in the stimulation field. In addition, the validity of the model was confirmed with experiments and simulations prior to experiments were carried out for the elimination of blind estimation [51]. Li et al. solved the PA Helmholtz equation by using the separation of variables method in spheroidal coordinates for a spheroidal droplet which is a more realistic geometric shape than a sphere [52]. When the PA detection of RBCs and their nuclei are considered, this derivation provides an opportunity to obtain a numerical simulation for practical uses of analyzing experimental data.
In this paper, we establish a PA wave expression for RBCs by using an extended Navier-Stokes equation when they travel into a microvessel. Our approach takes into consideration the interaction of cells to investigate RBC aggregation and its effect on blood rheology. We model the cells as spheres and the interaction forces between two cells cause the formation of aggregate clusters. We also analyze PA signal characteristics in the dynamic modes by means of RBC aggregation associated with complex hemorheological changes as well as viscosity. The motion of deformable cells can be described by the immersed boundary method based on Navier-Stokes equation. In order to evaluate the degree of aggregation which is related to their surface structure and their shape, the intercellular interaction is represented by the Morse type potential function [53]. Using a Green's function approach, an extended PA wave equation including a modified Morse type potential is solved analytically in detail. We obtain the PA pressure wave generated by normal and pathological RBCs associated with distinctive parameters such as a coefficient of surface energy and a scaling factor controlling the interaction decay behavior. The PA signal not only includes determinant parameters for the deformability of the cell such as plasma viscosity, density, and RBC aggregation but also contains laser parameters such as beamwidth, pulse duration, and repetition rate since the temporal profile is considered Gaussian rather than the Dirac delta function. Therefore, the approach given in this work can be useful to analyze the deformability of RBCs by taking advantage of the PA phenomena. The study of the mathematical modeling and the experimental design and detection of PA waves in micro-channels can be helpful for an appropriate choice of therapy.

Combined Navier-Stokes and heat transfer equations for a compressible flow
In this section, we present an expression for PA waves resulted from RBCs which travel into microvessels by solving the extended Navier-Stokes equation including RBC aggregation kinetics. The schematic of our model is shown in Fig. 1. The red blood cell is radiated with a pulsed laser. The PA waves resulted from the red blood cell are received by the detector. The interaction between the red blood cells is depicted by the arrows.
Absorption of light leads to a thermal expansion where ρ 0 is the mass density, c p is the constant-pressure heat capacity per unit mass, κ is the thermal conductivity, T(r, t) is the temperature rise at position r and time t, and H(r, t) is the heating function defined as the optical energy deposited per unit time and per unit volume [54,55]. As a result of the thermal expansion, acoustic waves are generated. In this process, the mass density obeys the continuity equation where ρ, v, and β are the mass density, the velocity field, and the coefficient of the volume thermal expansion, respectively [55,56]. The velocity field v holds for the Navier-Stokes equation where p is the pressure, η is the shear viscosity, and ζ is the bulk viscosity [57]. The thermal diffusion can be neglected since the pulse duration is very short so that the combination of Eq.
(1) with Eq. (2) gives Under the assumptions of small amplitude photoacoustic wave and sufficiently small source power, the pressure and density can be written as p = p 0 + δp and ρ = ρ 0 + δρ where p 0 and ρ 0 are the equilibrium pressure and density with δp p 0 and δρ ρ 0 , respectively [57]. After the linearizations of the Eqs. (3) and (4) under the assumption of small amplitude photoacoustic wave with constant p 0 and ρ 0 , we obtain and In order to simulate the blood flow in a microfluidic device for the deformation of a single red blood cell in a microvessel, the relationship between the mass density and the velocity field in the fluid can be obtained by solving the extended Navier-Stokes equation with a body force denoted by f(r, t) and where ν = (ζ + 1 3 η) is the blood viscosity. Here, F n is the aggregation force between two cells, n is the index of RBCs and N is the total number of interacting RBCs. The interaction force between RBCs, F n is modeled by the Morse potential, and ∆ is an interpolation function (interpolation kernel) which enables to model the interaction between the immersed structure and the fluid. Here, Y n is the position of the center of the red blood cell and r is the flow position. The cells aggregate slowly and an equilibrium configuration occurs when the intercellular interaction is balanced via the membrane forces after the deformation of the RBC membrane in microfluidic devices [58,59].
Combining Eqs. (5) and (7), the PA wave resulted from a collection of erythrocytes (which is considered as an optical absorber) can be obtained by solving the following equation where c s is the speed of sound, Γ = βc 2 s c p , δp = c 2 s δρ, and η = (ζ + 4 3 η). Thus, Eq. (9) can be written in the following form where S(r, t) is the source term which can be described by The first and second source terms on the right-hand side of Eq. (11) are due to the heat-producing, optical absorption within the two-phase sphere composed of one absorber (red blood cell) and fluid (blood plasma). Here, the fluid and the solid occupy a domain, but they do not intersect; thus, we model the fluid-structure interaction with the immersed boundary method and the boundary integral method. Moreover, the basic behaviour of the interaction forces between two RBCs is illustrated in Fig. 2. The third source term on the right-hand side of Eq. (11) is due to the interaction force between the two RBCs. Therefore, we can write the source term as S(r, t) = S 1 (r, t) + S 2 (r, t) + S 3 (r, t).

the wave equation can be expressed in the frequency domain
where The solution of Eq. (9) can be found by utilizing the Green's function approach [60]. The Green's function for Eq. (12) is given by where [60]. The heating function can be decomposed into the spatial and temporal parts as follows H(r, t) = A(r)I(t).
The heating due to the pulsed laser beam is very localized in both time and space so that the temporal and spatial parts of the heating function can be expressed by Gaussian functions. The radial part of the heating function can be written as Here, p 0 (r) is the initial pressure rise just after the laser pulse and it is defined by where c κ is the isothermal compressibility [61,62]. This initial pressure rise can be described by a Gaussian function as can be seen in the following section [62]. The temporal profile can be described by where τ is the pulse duration of the laser. Substituting Eqs. (15) and (17) into Eq. (14) leads to Applying the Fourier transform to S 1 yields Moreover, the second source term in Eq. (11) can be written as follows The Fourier transform of S 2 (r, t) is The function f(r, t) in Eq. (9) can be expressed by where n is the index of the particle and N is the number of absorbing particles (RBCs) referring to the radius of clusters of RBC aggregates. According to the immersed boundary method, we consider an compressible three-dimensional deformable structure in a RBC membrane immersed in an compressible fluid domain. In this regard, the combination of the fluid and structure is possible, provided that they do not intersect [59,[63][64][65]. Based on this method, the contribution of each RBCs whose center of mass is located at Y n to the flow at position r (which is considered as a distance between the source and the detector) is smoothed by a Gaussian distribution kernel. The choice of Gaussian kernel is a good assumption to study problems involving fluid-structure interactions in which an elastic structure is immersed in a viscous fluid. This method describes the fluid-membrane interaction between the flow field and deformable cells where X = (r − Y n (t)), and h is the standard deviation of the kernel, h = R √ π . Here, r varies between 0 and R and it is chosen as Y n .
The aggregating force between two cells represented by F n in Eq. (22) is modeled by a Morse type potential function. The Morse type potential function was also employed to describe the interactions of RBCs [53,66]. The parameters of the Morse potential function for intercellular interaction are summarized in Table 1 [67][68][69][70][71] and Table 2 [65,72] (In Table 1, the physical properties are mapped onto the dimensionless properties for the consistency of the simulation with the experiment. The scaling procedure relates the model's non-dimensional units to the physical units [70].). Figure 2 shows the aggregation forces corresponding to the various aggregation strengths which are summarized in Tables 3 and 4. These red blood cell aggregation forces calculated from the selected parameters show a good agreement with the experimental ones ( Fig.  2) [73][74][75]. These experimental forces are in the range of 2-12 pN [73], 14-23 pN, and 43-169 pN [74,75]. The measurements of interaction forces between red blood cells in aggregates and the measurements on dextran-induced aggregation of red blood cells were obtained by optical tweezers and atomic force microscopy-based single cell force spectroscopy, respectively.
We use the following Morse function to detect the acoustic pressure of non-aggregated and aggregated erythrocytes, where β, D, r 0 , and d are a scaling factor controlling the decay behavior of interaction of RBC, the coefficient of surface energy, the zero force separation, and the local distance between two surface elements of the cells, respectively. It is important to note that the number of absorbing particles (RBCs) depends on the radius of clusters of RBC aggregates (R). This number is expressed in terms of the radius of the defined region of interest since the amplitude of PA signal depends on radius of the absorber. In our case, N = 1. In our calculations, the interaction distance of the Morse potential regarding the immersed boundary method is calibrated with the distance between the source and detector in the case of the third source since our model couples the Navier-Stokes equation with the cell interaction and the bio-heat transfer equation. To see the contribution coming from the aggregation process on the photoacoustic wave, the interaction distance is taken as the distance between detector and the absorber. For attractive forces, the value of the scaling force is negative when r 0 < d. Conversely, if the value of the force is positive when r 0 > d, the force is indicated by a repulsive force. The attractive force is crucial to represent an intercellular interaction through the Morse-type potential energy function. Under this intercellular interaction, the aggregation occurs either by increasing D or by decreasing β, which is shown in Figs. 2-4. Moreover, zero force length may raise the depletion thickness as well as equilibrium distance of the RBCs in aggregates. In

Solution of the extended photoacoustic wave equation for a Gaussian radial absorption profile
In this section, we describe the radial profile by a Gaussian function with the standard deviation or beamwidth of the laser σ where θ is the Heaviside step function. Here, the initial pressure p 0 is created inside the object so that the initial pressure distribution can be defined as Eq. (25). The following integral can be written as a summation of three terms corresponding to three different sources as follows where S=S 1 +S 2 +S 3 . The contribution coming from the first source is where | r − r |= r 2 + r 2 − 2rr µ , µ = cos θ , and So the solution due to the first source becomes Using the Taylor-series expansion we evaluate the integral in Eq. (28) [60, 76] The exponential term can be further expanded [60,76,77] exp iω(r ± r ) Thus, we obtain [76,77] Here, it is important to note that after a series of numerical calculation, we see that the contribution coming from the high order terms is negligible. This case can also be analyzed analytically.
For large values of ω, the compactly supported function exp(−τ 2 ω 2 /2) dominates the higher order terms resulted from the series expansion. On the other hand, the physical parameters pulse duration (τ) and B = η /(ρ 0 c 2 s ) are also on the orders of 10 −9 and 10 −12 , respectively so that the value of the integral is already very small for small values of ω. Therefore, based on the relatively small values of τ and B and the predominance of exp(−τ 2 ω 2 /2) term for large values of ω, the influence of the integral due to the higher order terms is negligible. So that the p 1 (r, t) can be simplified as follows where and Substituting the results of the integrals J 1a and J 1b into Eq. (33) enables to find the pressure due to S 1 (r, t) where erf(x) is the error function. It is important to note that the value of the integral J 1b is sufficiently small compared to the integral J 1a .
The solution of the PA equation due to the second source in time domain can be evaluated by using the Taylor series expansion in the following p 2 (r, t) expression similar to the calculation of p 1 (r, t) Hence, p 2 simplifies to The first integral with respect to ω in Eq. (38) can be calculated by using the residue theorem. Applying the residue theorem at z=0 The contribution coming from integral J is zero so that p 2 (r, t) becomes We can express p 2 (r, t) in the following form and Moreover, the third PA wave, p 3 (r, t) due to the third source corresponding to aggregation force is The first integral over ω in Eq. (45) has the same form as shown in Eq. (38), and it can be evaluated by the contour integration method. We can rewrite Eq. (45) in the following form which yields the pressure wave of non-aggregated and aggregated erythrocytes where and By evaluating the integral in Eq. (47), we obtain Here, the Heaviside step function, θ, is introduced to take into consideration the propagation time when the observation point is outside the PA absorber [61,62]. The integrals J 3a and J 3b contribute significantly to the PA signal. We note that the calculation of the integral J 3a is very similar to the one of the integral J 3b .

Numerical/simulation parameters
The behavior of a normalized PA wave produced by a RBC as a function of normalized time is illustrated in Fig. 5 for two different pulse durations (τ=1 and 5 ns) and beamwidths (σ = 6 and 8 µm). The PA wave stemming from the first source is calculated from Eq. (33). In Fig. 5, the effects of the laser beamwidth and pulse duration on the PA signal are observed. The density, the speed of sound within the cell, and the radius of absorber are taken as ρ = 1000 kg/m 3 , c s =1520 m/s, R=8 µm [78], respectively. Fig. 6 shows the power spectrum of the PA signal given in Fig. 5. Even though the change in the PA amplitude is small in the time domain, it becomes more visible in the frequency domain. Fig. 7 illustrates the effects of the laser beamwidth and pulse duration on the PA signal for the detector position (2 mm away from the center of the absorber [79,80]) regarding the feasibility of experimental realization in the presence of the first source. Fig. 8 and Fig. 9 show the effects of the plasma viscosity on the normalized PA wave, where an ultrasonic detector is located at r =2R and 500 µm, respectively [79,81]. In Figs. 8-15, the shear and bulk viscosities are taken as η = 1.2, 1.3 mPa.s, ζ = 4.5, 5, 6, 18.4, and 22.9 mPa.s, respectively [74,[82][83][84] as can be seen in the following section.

Results
The integrals J 1a , J 1b , J 2a , and J 2b given by Eqs. (34), (35), (43), and (44) are evaluated numerically and the pressure wave is found as a summation of p 1 and p 2 corresponding to the first and second sources, respectively. However, the contribution coming from the second source is negligible compared to the first one. The pressure waves from the first and second sources can be applied for the suspension of RBCs in an compressible fluid. Our results show that the magnitude of PA signal does not change with different plasma viscosity values as shown in Fig.  8, as expected, in the presence of first and second source terms. Moreover, Fig. 9 shows that the amplitude of PA wave is decreasing with the distance from absorber, r=500 µm when fiber optic sensing systems are considered [79]. The aggregation interaction is described by the Morse potential consisting of a short-range repulsive force when d < r 0 and a long-range attractive force for d > r 0 [66,85] The corresponding force, F = −∂E/∂r, is obtained from the Morse potential. In our calculations, the energy depth D is used in units of Joule and given in Tables 2 [65,72].   In this paper, the force is approximated with the several parameters involved in the detection of the PA signal, and the immersed boundary method is used to eliminate the effect of the fluid-cell interaction.
Different aggregating conditions are created by changing the surface energy (D), the scaling factor which is related to the thickness of the depletion layer or interactive distance (β), and the zero-force length (r 0 ) which are shown in Tables 3 and 4. In these conditions, the aggregation strength is directly related to amount of surface energy, which is expressed by the interaction force between RBCs. The parameters of the Morse potential in Table 3 used in this study were taken from the works presented by Yazdani et al. [65] and Fenech et al. [72], and the ranges of forces shown in Fig. 2(a) were fitted with the experimental data [73]. Moreover, the ranges of forces in Table IV were taken from studies presented by Neu et al. [75] and Steffen et al. [74], as shown in Fig. 2(b). In Fig. 10, the PA pressure wave generated by three sources increases with RBC aggregation (case 1 represents the high aggregation among three cases in Table 3, which is shown in Fig. 2). Moreover, our results demonstrate that high whole-blood viscosity has a role to play in erythrocyte aggregation since the amplitude of PA signal increases when the whole-blood viscosity is higher, as shown in Fig. 10 for case 1. When we compare the three cases, cases 2 and 3 show a lower aggregation. The influences of laser parameters and the levels of RBC aggregation on PA wave are shown in Figs. 11 and 12. These Figs. 11 and 12 also show that the short pulse duration results in large wave amplitude, and the amplitude decreases as beamwidth diminishes for these cases. The aggregation for cases 2 and 3 is low and the force in the case 2 is close to the one in case 3. We examine the cases of Dextran 70 and Dextran 150 in Figs. 13 and 14. In Fig. 15, we also analyze the cases 1 and 3 allowing to differentiate and characterize their aggregation levels by observing the PA signal amplitude regarding the chosen force corresponding to the values of function of separation between cell surfaces (β(d − r 0 )) as shown in Fig. 13. Here, the same points (β(d − r 0 )) for each case correspond to the various aggregation forces as shown in Fig. 14(d). Our theoretical simulations show that PA signal amplitude increases linearly with increasing aggregation force, which are in excellent agreement with literature data [15,50,88,89]. We use the measured interaction forces for types of Dextran solutions inducing the aggregation in order to validate our analytical approach. Moreover, our results suggest that PA can be used to differentiate red cell aggregation process in detail by using a known interaction force between cells.

Discussion
The PA imaging may have a great potential for anatomical, functional and molecular imaging of RBCs. In this work, the level of aggregation of RBCs in human blood is assessed and the morphology of blood cells is studied for the detection of RBCs using PA waves obtained analytically. To see the feasibility of the proposed extended photoacoustic transport model including a modified Morse type potential function, an analysis of the amplitude of PA signal is conducted by using the measured aggregation forces [73][74][75]. We compare our theoretical results with the previous experimental and simulation studies for the PA detection of RBC aggregation and see that our results are in good agreement with the literature [46, 47, [73][74][75]. Our numerical results show that the PA signal amplitude increases with the level of aggregation. Some experimental-simulation studies presented that the PA signals could be used to identify the presence of varying degrees of aggregation [15,50]. However, many other factors such as effects of flow, intervening tissues, ultrasound transducer receiver bandwidth are essential for the experimental confirmation of this study [50]. In order to analyze PA radiofrequency (RF) signals produced by RBC aggregates, Hysi et al. reported the feasibility of detecting these PA RF spectral changes by incorporating a finite transducer in detection of PA signals in his theoretical and experimental study [15]. In our model, the combined Navier-Stokes and PA equation system is solved to find the PA wave by utilizing the immersed boundary method for fluid-structure interaction. Furthermore, our analytical approach shows the effect of the whole blood viscosity on the aggregation processes in the case of third source due to the intercellular interaction force and it also allows to analyze each source separately (If we want to include the interaction force for more than one RBC, then the heating function (H) would be expressed as the summation of the contributions resulted from each RBC.).
In this paper, we solve a combined heat conduction and Navier-Stokes equation system utilizing the Morse-type potential energy function for the characterization of human red cell. The fluid flow containing RBCs in microvessel is modeled by an extended Navier-Stokes equation with the immersed boundary method proposed by Peskin et al. [90]. The immersed boundary method has been utilized to model the flow of RBCs in blood plasma or to describe the blood plasma and RBC membrane interaction between the flow field and deformable cells [59]. We modified the immersed boundary method with respect to the location of PA detector and radius of PA source in order to obtain the PA wave coming from RBC aggregation. Furthermore, our solution makes it possible to observe the effect of laser parameters on the PA wave because the temporal profile is taken as Gaussian rather than a point like Dirac delta function. To validate our method with another one presented by Erkol et al.'s method [62], we also show that the PA signal decreases when the beamwidth decreases and shorter pulse duration leads to a larger wave amplitude as expected. In this regard, adjustment of laser parameters can be very useful for the application of the detection of RBCs.
Another advantage of our method is that it provides comprehensive solutions since the spatiotemporal profile of the laser is taken Gaussian (rather than point source) and RBC aggregation and plasma parameters are taken into consideration. It is well known that blood rheology and viscosity are related to hematocrit (Ht), RBC aggregation, shear induced deformation of RBCs, and plasma viscosity [91][92][93][94]. In recent years, several studies have demonstrated that the level of RBC aggregation increases when the plasma fibrinogen concentration increases, which shows the capability of PA for dynamic monitoring of rheological parameters in circulating blood [71,[95][96][97]. Our results show that aggregation is associated with the whole blood plasma viscosity after RBC aggregation because the pressure wave amplitude is promoted by the alteration of whole blood viscosity in the case of the third source due to the intercellular interaction force. The plasma viscosity may increase up to five-fold in an unfavorable situation ranging from inflammatory diseases to plasma cell dyscrasias [98,99]. Therefore, the plasma viscosity is significant in specific cases only although this influence is relatively small in general cases [100][101][102]. For an application, we studied the red blood cell aggregation in a microcirculatory system. Our detailed analytical solution of the extended photoacoustic transport model combined with a modified Morse type potential function accounts for the dynamics of aggregate formation. Hence, the distinctive features of the photoacoustic waves enable to understand diseased or damaged RBCs having abnormal shapes and sizes, which can result from a variety of diseases. Our detailed solution of PA wave to detect RBC aggregation is obtained by solving the combined Navier-Stokes and bio-heat transfer equations. This study also presents the spectral features of the PA signals by considering the interaction between RBCs modeled by the Morse potential. We think that our work can be helpful to improve hybrid simulation approach for not only the characterization of red blood cell morphology but also for cancer cell lines in the microchannel flow. To describe and analyze the mechanical processes inside a microfluidic device, our hybrid model can pave the way toward the development of new physical models in the rapidly advancing field of lab-on-chip technology and photoacoustics. As a future work, in order to analyze the influence of the deformation of the membrane on pressure wave, we will implement the spring network approach beside the Morse potential.