Transport dissipative particle dynamics model for mesoscopic advection-diffusion-reaction problems

Many biological processes take place at the cellular and subcellular levels, where the continuum deterministic description is no longer valid and hence stochastic effects have to be considered.¹ To this end, mesoscopic methods with stochastic terms are attracting increasing attention as a promising approach for tackling challenging problems in bioengineering and biotechnology.² As one of the currently most popular mesoscopic methods, dissipative particle dynamics (DPD) drastically simplifies the atomistic dynamics by using a single coarse-grained (CG) particle to represent an entire cluster of molecules,³ and the effects of unsolved degrees of freedom are approximated by stochastic dynamics.⁴ Similarly to the molecular dynamics (MD), a DPD system consists of many interacting particles and their dynamics are computed by time integration of Newton’s equation of motion.⁵ However, in contrast to MD, DPD has soft interaction potentials allowing for larger integration time steps. With larger spatial and temporal scales, DPD modeling can be used to investigate hydrodynamics in larger systems, which are beyond the capability of conventional atomistic simulations.⁶ 

DPD was initially proposed by Hoogerbrugge and Koelman³ for simulating the mesoscopic hydrodynamic behavior of complex fluids. The interactions between DPD particles occur pairwise so that the total momentum of the DPD system is strictly conserved. Moreover, since these interactions depend only on the relative positions and velocities, the resulting DPD fluids are Galilean invariant.⁵ By using the Fokker-Planck equation and applying the Mori projection operator, Español⁷ and Marsh et al.⁸ reported that the hydrodynamic equations of a DPD system recover the continuity and Navier-Stokes equations. In other words, DPD can be considered as a particle-based Lagrangian representation of the continuity and Navier-Stokes equations at mesoscopic level. Therefore, DPD can provide the correct hydrodynamic behavior of fluids at the mesoscale, and it has been successfully applied to study DNA suspensions,^9,10 blood flows and cell interactions,^11,12 platelet aggregation,¹³ and nanoparticles across lipid membranes.¹⁴ 

However, the classic DPD system has only equations governing the evolution of density and velocity fields but no evolution equations for describing the concentration fields, which precludes DPD from modeling problems involving diffusion-reaction processes, i.e., two of the most fundamental processes in biological systems.¹ Proteins in an aqueous solution diffuse in a living cell due to Brownian motion, and some collisions of appropriate proteins may lead to chemical reactions. Moreover, most biological processes depend on the concentrations of specific proteins, ions, or other biochemical factors.¹⁵ For example, blood coagulation, a crucial process involved in thrombus formation and growth, can be modeled by advection-diffusion-reaction (ADR) processes of 25 biological reactants.¹⁵ As another example, the sickle-cell disease is initiated by the deoxygenation-induced polymerization of sickle hemoglobin in the red blood cells, and the polymerization dynamics relies heavily on the concentration of oxygen.¹⁶ Therefore, it is important to include the ADR equation in the DPD model when diffusion and reactions are present.

Since the ADR equation has been extensively investigated by continuum solvers, in which smoothed particle hydrodynamics (SPH) is a Lagrangian discretization of the continuum equations,¹⁷ it is necessary to clarify the differences of the mesoscopic stochastic method, DPD, from continuum deterministic methods, i.e., SPH. Although the hydrodynamic behavior of a DPD fluid can recover the Navier-Stokes equation, the DPD model does not come from continuum equations. Instead, DPD is usually considered as a CG MD model, and its formulations can be obtained from applying the Mori-Zwanzig projection to MD systems.⁴ Table I presents a qualitative comparison between DPD and SPH methods. It shows that the DPD method has advantages on modeling complex fluids and materials, which may even not have a known constitutive equation.

In the present paper, we develop a transport dissipative particle dynamics (tDPD) model for mesoscopic problems involving ADR processes. Similarly to the classic DPD method, each tDPD particle is associated with extra variables for carrying concentrations in addition to other quantities such as position and momentum. Using the same strategy for discretization of the heat conduction equation,^18,19 the transport of concentration can be modeled by a Fickian flux and a random flux between particles.^20,21 An analytical formula is proposed to relate the tDPD parameters to the effective diffusion coefficient. Furthermore, we develop a methodology for implementation of the correct Dirichlet and Neumann boundary conditions in tDPD simulations. It is worth noting that the particle-based tDPD method satisfies the conservation of concentration automatically and provides an economical way to solve ADR equations with a large number of species. For instance, in a problem with 23 ADR equations involving 25 chemical species, the additional computational cost of tDPD is the same as the flow solver, unlike the continuum model requiring 23 Poisson/Helmholtz solvers making the additional computational cost over ten times higher than a pure Navier-Stokes solver.²² 

The remainder of this paper is organized as follows. In Section II, we describe the details of tDPD formulation and its parameters, including an analytical formula to relate the mesoscopic concentration friction to the effective diffusion coefficient (Section II A) and a methodology for imposing Dirichlet and Neumann boundary conditions (Section II B). In Section III, we validate the tDPD model with one-dimensional diffusion (Section III A) and two-dimensional advection-diffusion (Section III B) and present an application of the tDPD model to ADR processes involving 25 species in coagulation and fibrinolysis in flowing blood (Section III C). Finally, we conclude with a brief summary in Section IV.

tDPD is an extension of the classic DPD framework. With extra variables for describing concentration fields, the time evolution of a tDPD particle i with unit mass is governed by the conservation of momentum and concentration, which is described by the following set of equations:

where t, r_i, v_i, and F_i denote time, and position, velocity, force vectors, respectively. $F i ext$ is the force on particle i from an external force field. The pairwise interaction between tDPD particles i and j consists of the conservative force $F ij C$⁠, dissipative force $F ij D$⁠, and random force $F ij R$⁠, which are expressed as⁵ 

where r_ij = |r_ij| with r_ij = r_i − r_j, e_ij = r_ij/r_ij is the unit vector from particle j to i, and v_ij = v_i − v_j the velocity difference. ω_C(r), ω_D(r), and ω_R(r) are the weight functions of F^C, F^D, and F^R, respectively. Also, a_ij is the conservative force parameter, γ_ij the dissipative coefficient, and σ_ij the strength of random force. Moreover, C_i represents the concentration of one species defined as the number of a chemical species carried by a tDPD particle i and Q_i the corresponding concentration flux. Since tDPD particles have unit mass, this definition of concentration is equivalent to the concentration in terms of chemical species per unit mass. Then, the volume concentration, i.e., chemical species per unit volume, is ρC_i where ρ is the number density of tDPD particles. We note that C_i can be a vector C_i containing N components, i.e., {C₁, C₂, …, C_N}_i when N chemical species are considered.

Based on Fick’s law,²³ the diffusion driving force of each species is proportional to the concentration gradient, which corresponds to a concentration difference between two neighboring tDPD particles. Therefore, in the tDPD model, the total concentration flux on particle i accounts for the Fickian flux $Q ij D$ and random flux $Q ij R$⁠, which are given by

where κ_ij and ϵ_ij determine the strength of the Fickian and random fluxes. The symbols ξ_ij in Eq. (5) and ζ_ij in Eq. (7) represent symmetric Gaussian random variables with zero mean and unit variance.⁵ A similar DPD transport model was proposed in Ref. 20 by drawing an analogy with the DPD model of heat conduction.¹⁸ It is worth noting that the parameter κ plays the same role of dissipating the concentration differences between tDPD particles as the dissipative coefficient γ does on momentum. To avoid confusion with γ, we use “Fickian friction coefficient” for κ to denote the dissipative effect of concentration. In general, the Fickian friction coefficient κ is an N × N matrix when the interdiffusivities of N different chemical species are involved. In simple cases, i.e., considering N chemical species in dilute solutions and neglecting the interdiffusivities of different species, κ becomes a diagonal matrix.²⁴ 

The random force $F ij R$ and the random flux $Q ij R$ correspond to the irreversible dissipative force $F ij D$ and dissipative Fickian flux $Q ij D$⁠, respectively. The fluctuation-dissipation theorem (FDT) is employed to relate the random terms to the dissipative terms^25,26

where k_B is the Boltzmann constant, T the temperature, and m_s the mass of a single solute molecule. Also, C_i and C_j are the concentrations on particles i and j, respectively. Simplified derivations for obtaining Eq. (9) are presented in Appendix  A. Since m_s is small with respect to the mass of a tDPD particle m, which is usually chosen as mass unit, the magnitude of ϵ is small. Therefore, the contribution of the random flux $Q ij R$ to the total diffusion coefficient D is negligible unless m_s becomes comparable to m in nanoscale systems. In the present work, we consider m_s/m ≪ 1 and neglect the contribution of the random flux $Q ij R$ to the total diffusion coefficient.

The third part in Eq. (2) $Q i S$ represents a source term of concentration due to chemical reactions. Since the source term $Q i S$ relies on specific chemical reaction equations, here we do not present a general form of $Q i S$⁠. In the applications of the tDPD model in Section III, the reaction equations of 23 chemical species are involved and the expressions of $Q i S$ will be presented. We note that a random rate of chemical reactions corresponding to the deterministic reaction rate $Q i S$ can be also considered by applying FDT to particular chemical reaction equations.²⁶ Although we do not consider the random reaction term in the present work, it can be readily included in a fluid system of interest involving specific chemical reactions.

The macroscopic properties including viscosity and diffusivity of a tDPD system are output properties rather than input parameters. Since the stochastic forces on tDPD particles yield random movements, the effective diffusion coefficient D is the result of contributions from the random diffusion D^ξ and the Fickian diffusion D^F. In general, the random contribution D^ξ is a combined result of the random movements of tDPD particles and random flux $Q ij R$ in Eq. (7). However, the variance of random flux $Q ij R$ has a small prefactor $m s 2$ as given by Eq. (9). The contribution of the random flux $Q ij R$ to D^ξ is negligible, which has been confirmed by Kordilla et al.²¹ Thus, in the mathematical derivations in this section, we consider that D^ξ is due to the random movements of tDPD particles.

Considering a tDPD system in which the particles are governed by the interactions given by Eqs. (3)-(5), the diffusion coefficient D^ξ induced by the random movements of tDPD particles can be calculated by⁵ 

where r_c is the cutoff radius for forces. Also, the macroscopic diffusion coefficient D^F due to Fickian concentration flux can be computed by²⁷ 

where κ is the Fickian friction coefficient and r_cc is the cutoff radius for concentration flux. The forces between tDPD particles include indirect effects of interactions while the exchanges of chemical species need direct contact between particles. Therefore, we suggest to use a cutoff radius r_cc less than or equal to r_c. Specifically, when the radial distribution function of ideal gas g(r) = 1.0 is employed, both D^ξ and D^F can be evaluated analytically and we have

where s₁ and s₂ are the exponents of the weighting functions given as ω_D = (1 − r/r_c)^s₁ and ω_DC = (1 − r/r_cc)^s₂, respectively. Equation (12) provides a relationship between the macroscopic effective diffusion coefficient D (which can be experimentally measured) and parameters in the tDPD model. Specifically, Eq. (12) defines scaling of D^ξ and D^F with temperature and the particle size (i.e., r_c and r_cc).

Equation (12) indicates that the effective diffusion coefficient D is a linear function of the parameter κ, and the minimum value of the effective diffusion coefficient D_min = D^ξ is obtained at κ = 0. We note that tDPD simulations should have smaller time steps for larger κ to avoid instability issues. Considering a tDPD particle i with its concentration C_i surrounded by other particles with concentration of C_i + ΔC, the concentration flux on particle i can be computed by a integral of $Q ij D$ in a spherical domain with radius of r_cc. Then, the variation of concentration on particle i during one time step in an explicit time integration must be smaller than ΔC, which results in a necessary condition for limiting the time step Δt given by

By assuming the radial distribution function g(r) = 1 corresponding to ideal gases,^5,19 Eq. (12) can be used to predict the effective diffusion coefficient of a tDPD system roughly, but the accurate value of D can only be obtained by computations in tDPD simulations. For a tDPD fluid with constant diffusion coefficient, the ADR equation is given by

where D is the diffusion coefficient and Q^S a source term. For steady state problems, Eq. (14) is simplified to D∇²C = − Q^S, which is the same as the governing equation of Poiseuille flow driven by a body force υ∇²V = − g. Therefore, we can use a method analogous to the reverse Poiseuille flow²⁸ for computation of the effective diffusion coefficient.

To obtain the accurate value of the effective diffusion coefficient, we perform tDPD simulations in a computational domain of 30 × 6 × 30 with periodic boundary conditions. The parameters in the tDPD system are defined as ρ = 4.0, k_BT = 1.0, a = 75k_BT/ρ, γ = 4.5, σ² = 2k_BTγ, r_c = r_cc = 1.58, ω_C(r) = (1 − r/r_c), $ω D ( r ) = ω R 2 ( r ) = ( 1 − r / r c ) 0 . 41$⁠, and ω_DC(r) = (1 − r/r_cc)² and κ ranges from 0 to 10. As shown in Fig. 1(a), the fluid system is subdivided into two equal domains in z-direction, and a concentration source +Q^S is applied in the domain of z > 0 while a concentration sink with same magnitude −Q^S is applied in the other domain z < 0. Because of the periodic boundary conditions, the concentration of the tDPD fluid is constrained to be invariable at the plane z = 0. When the diffusion coefficient D is constant, the steady state solution of the concentration profile is given by

where Q^S = 1.0 × 10⁻³ is the source term, C₀ = 1.0 the initial concentration of the tDPD system, and d = 15.0 the half length of the computational domain in z-direction. The concentration profiles shown in Fig. 1(a) were obtained by subdividing the computational domain into 51 bins along the z-direction. The tDPD simulations were run for 500 time units and the symbols in Fig. 1(a) are the concentration profiles of the last time step in the simulations. Then, the effective diffusion coefficient can be determined by fitting the concentration profile with the analytical solution given by Eq. (15). The solid lines in Fig. 1(a) are the best-fitting parabolas. It is obvious that the effective diffusion coefficient D can be significantly changed by varying the Fickian friction coefficient κ.

Figure 1(b) shows the dependence of the effective diffusion coefficient D on the Fickian friction coefficient κ. As indicated by Eq. (12), D should be a linear function of κ, which is verified by the computed results shown in Fig. 1(b). Inserting the tDPD parameters in the approximate theoretical prediction, Eq. (12) gives D = 0.0195 + 0.786κ, while the linear fitting of computed data shown in Fig. 1(b) gives D = 0.0187 + 0.771κ. Hence, the analytical formula of Eq. (12) predicts roughly the diffusion coefficients with approximately a 5% deviation from the computed results. In practice, the accurate value of D should be computed by performing tDPD simulations with the method analogous to the reverse Poiseuille flow as shown in Fig. 1(a).

Boundary conditions are crucial for the investigation of diffusion-reaction processes in wall-bounded systems. Usually, defining stationary particles to represent solid objects is a common treatment in classic DPD simulations.²⁹ However, those solid walls made up by discrete frozen particles induce unwanted temperature and density fluctuations in the vicinity of the walls.^30,31 Alternatively, Lei et al.³² reported that the effects from a solid wall can be evaluated by integrating the interactions between fluid and solid particles, and hence, the presence of solid wall boundaries was replaced by effective boundary forces. Subsequently, Li et al.¹⁹ extended this method to non-isothermal systems for applying thermal boundary conditions and they found negligible temperature and density fluctuations near the wall boundaries. In the present work, we adopt the idea of effective boundary fluxes to impose Dirichlet and Neumann boundary conditions for concentration in the tDPD systems.

Since the fluid particles do not penetrate into wall boundaries, the random movements of fluid particles do not have any contributions to the boundary concentration flux. Therefore, the effective boundary concentration flux is induced only by the Fickian flux. For a fluid particle i in the vicinity of a wall boundary, the effective Fickian concentration flux on particle i from the wall can be calculated¹⁹ by

where C_i is the concentration of particle i, C_w the expected concentration at the wall surface, and h the distance of the particle i away from the wall surface. Here, φ(h) is a function of h, which is defined as

Equation (16) reveals that the boundary concentration flux is determined by the concentration difference between particle i and the wall, and also their distance.

For a given tDPD system with parameters ρ = 4.0, k_BT = 1.0, a = 18.75, γ = 4.5, σ = 3.0, r_c = 1.58, ω_C(r) = (1 − r/r_c), and $ω D ( r ) = ω R 2 ( r ) = ( 1 − r / r c ) 0 . 41$⁠, the radial distribution function g(r) of the tDPD particles is displayed in Fig. 2(a). Then, the function φ(h) can be evaluated through Eq. (17) with ω_DC(r) = (1 − r/r_cc)² and r_cc = r_c. We note that though r_cc is usually set to be r_c for convenience, r_cc can be different from r_c. The dependence of φ(h) on the distance h is presented in Fig. 2(b), which is approximated by φ(h) = 4.9562 × 10⁻²/h ⋅ [1 + 4.1517h/r_cc + 7.1697(h/r_cc)²](1 − h/r_cc)^4.09 for h < r_cc. As the distance h approaches to zero, the magnitude of φ(h) goes to infinity. Here, we use a truncation of φ(h) at small distances and set φ(h < 0.01r_cc) = φ(0.01r_cc).

To consider the effective flux along the normal direction of wall surface, we integrate the effect of concentration flux from the wall boundary and define a distance dependent function given by

The normalized Φ(h) is defined as $Φ ¯ ( h ) =Φ ( h ) / ∫ 0 r c c Φ ( x ) dx$⁠. Then, the integral of $Φ ¯ ( h )$ is equal to one. Using the computed g(r) shown in Fig. 2(a) and ω_DC(r) = (1 − r/r_cc)², the function $Φ ¯ ( h )$ can be obtained through Eq. (18), as shown in Fig. 2(b). The computed $Φ ¯ ( h )$ is approximated by $Φ ¯ ( h ) = ∑ i = 0 6 p i ( 1 − h / r c c ) 7 − i$ with the coefficients p = {22.745 45, − 62.6256, 55.1365, − 19.2584, 6.2662, − 0.502 42, 0.0101} for h ≤ r_cc and $Φ ¯ ( h ) =0$ for h > r_cc.

To impose a Neumann boundary condition dC/dn = Λ at a wall boundary, it is equivalent to apply a concentration flux Q^W = DΛ across the boundary. In practice, the flux Q^W is distributed onto the fluid particles in the vicinity of the wall weighted by $Φ ¯ ( h )$⁠. Let ρ be the number density of the fluid particles, the volume concentration is ρC because the concentration C in tDPD is defined as the number of a chemical species per particle. Then, any fluid particle i close to the wall obtains a concentration flux from the wall boundary given by

In addition to effective boundary concentration flux, the interactions between wall and fluid particles are also integrated and replaced by effective boundary forces. In particular, by using the boundary method proposed in Refs. 19 and 32, the effective conservative and dissipative forces are given by

where the unit vector e_n represents the normal direction of wall surface. With the computed g(r) as shown in Fig. 2(a), we have $f C ( h ) = ∑ i = 0 6 p i ( 1 − h / r c ) 7 − i$ with the coefficients p = {525.6666, − 1457.0412, 1485.7229, − 927.0881, 453.3782, − 7.2616, 0.177 58},

where v_ne_n = (Δv ⋅ e_n)e_n is the projected component of velocity difference Δv = v_i − v_wall onto the normal direction of wall surface, and v_τe_τ = Δv − v_ne_n is the tangential component of velocity difference Δv. Here, the symbols τ and n represent the tangential and normal directions of wall surface, respectively. Using the computed g(r) and ω_D(r) = (1 − r/r_c)^0.41, both γ_τ(h) and γ_n(h) can be computed through Eq. (21). In practice, the computed γ_τ(h) and γ_n(h) are approximated by γ_τ(h) = 10.3101/h ⋅ [1 + 3.6212h/r_c + 2.2090(h/r_c)²](1 − h/r_c)^3.45 and γ_n(h) = 20.7789/h ⋅ [1 + 2.1245h/r_c + 5.4566(h/r_c)²](1 − h/r_c)^2.30 for h < r_c, respectively.

In addition, the effective boundary forces include also the fluctuating part to satisfy the fluctuation-dissipation theorem, which is given by

where ξ_τ and ξ_n are Gaussian random variables with zero mean and unit variance. The variances of random forces are related to the dissipative forces by the fluctuation-dissipation theorem in the form of $σ τ 2 ( h ) =2 k B T γ τ ( h )$ and $σ n 2 ( h ) =2 k B T γ n ( h )$⁠. When the distance h approaches to zero, both γ_τ(h) and γ_n(h) become infinity. In practice, we set a maximum magnitude for γ_τ(h) and γ_n(h) given as γ_τ(0.05r_c) and γ_n(0.1r_c) to avoid instability of tDPD simulations induced by some extremely large random numbers.

The boundary method based on the effective boundary forces has been well validated by Lei et al.³² and Li et al.¹⁹ They have shown that the density and temperature fluctuations in the vicinity of wall surface are successfully eliminated by using those effective boundary forces. In the present paper, we will not recheck the density and temperature fields and will focus on the concentration fields instead.

The first test case is solving the diffusion equation with the initial condition given by the Heaviside step function. Considering the diffusion in a system containing two components in dilute solution, by neglecting the interdiffusivities of different species, the system is reduced to two uncoupled diffusion equations with diffusivities independent between species. Then, the diffusion equation takes the relatively simple form of the linear second-order partial differential equation

where C₁ and C₂ are the concentration of species 1 and 2, respectively. D₁₁ and D₂₂ are the corresponding diffusion coefficients. In general, the diffusion coefficient can be a function of concentration, time, and direction to consider heterogeneous diffusions. However, here we consider only the simplest case in which D₁₁ and D₂₂ are constant.

As shown in Fig. 3(a), it is a one-dimensional diffusion along x-direction. In practice, a tDPD system is defined in a computational domain of 100 × 10 × 10 with periodic boundary conditions in y- and z-directions and solid walls in x-direction. The concentration of the tDPD system is initialized as C₁ = 1.0 and C₂ = 0.0 in the domain of x < 0 and C₁ = 0.0 and C₂ = 1.0 in the domain of x ≥ 0. The solid walls provide no-slip boundary condition for velocity and no-flux boundary condition for concentrations. At wall boundaries, we apply the effective boundary forces given by Eqs. (20)-(22) and no concentration flux Λ = 0 shown in Eq. (19). The parameters of the tDPD system are given as ρ = 4.0, k_BT = 1.0, a = 75k_BT/ρ, γ = 4.5, σ² = 2k_BTγ, r_c = r_cc = 1.58, ω_C(r) = (1 − r/r_c), $ω D ( r ) = ω R 2 ( r ) = ( 1 − r / r c ) 0 . 41$⁠, and ω_DC(r) = (1 − r/r_cc)². Unless otherwise specified, all the simulations in the present paper are based on this set of parameters. Also, κ₁₁ = 0.04 results in D₁₁ = 0.05 and κ₂₂ = 0.62 yields D₂₂ = 0.50. The velocity-Verlet algorithm⁵ is utilized for the numerical integration of the tDPD equations of Eqs. (1) and (2). The time step Δt = 0.01 is used for Section III A and Δt = 0.001 for Sections III B and III C, which satisfy the necessary condition given by Eq. (13) for limiting the time step Δt.

For one-dimensional diffusion with initial condition of step function, the time-evolution of concentration profiles before concentrations diffuse to the walls on either side can be described by

where k = 1, 2 represent the two different species, and $erf ( x ) =2/ π ⋅ ∫ 0 x e − τ 2 dτ$ is the error function.

The concentration profiles shown in Fig. 3(b) were obtained by subdividing the computational domain into 100 bins along the x-direction. The local concentration is obtained by averaging all the particle information in each bin based on a single time step. The lines in Fig. 3(b) represent the theoretical solutions given by Eq. (24). The transient concentration profiles in the tDPD simulation are in excellent agreement with the theoretical solution, which validate the present tDPD model for solving the diffusion equation.

The second test case is used to validate the boundary method for implementation of the Dirichlet boundary condition. Figure 4(a) illustrates the initial condition C(x, 0) = 0 and boundary conditions. Similarly to the first test case, the tDPD simulation is performed in a computational domain of 100 × 10 × 10 with periodic boundary conditions in y- and z-directions and no-slip solid walls in x-direction. By solving a one-dimensional diffusion equation with boundary conditions of C(0, t) = 0 and C(100, t) = C₀, an analytical solution for the transient concentration profile can be obtained,

where β_n = nπ/L_x with L_x being the length of the computational domain in the x-direction, D the diffusion coefficient. The Fickian friction coefficient is κ = 5.0 and the corresponding diffusion coefficient is D = 3.87. The computational domain is divided into 100 bins along the x-direction for obtaining local concentrations. Figure 4(b) shows a comparison between the concentration profiles obtained using tDPD and theoretical solution Eq. (25) at several times including the steady state solution. The results are in good agreement, which validate the boundary method for imposing the correct Dirichlet boundary condition in the tDPD simulation.

The last one-dimensional test case is performed to validate the boundary method for implementation of the Neumann boundary condition. As shown in Fig. 5(a), it has a similar setup as the second test case but different wall boundary conditions. Considering a concentration flux at the left wall x = 0, we apply the Neumann boundary condition dC/dn = Λ = 0.01 at x = 0. Also, the right wall at x = 100 has a fixed concentration C(100, t) = 0. By solving a one-dimensional diffusion equation with boundary conditions of dC(0, t)/dx = Λ and C(L_x, t) = 0, we have the theoretical solution for the transient concentration profile given by

where β_n = nπ/(2L_x) with L_x being the length of the computational domain in the x-direction, D the diffusion coefficient. tDPD is used to simulate the time evolution of the concentration profile for D = 3.87. The Neumann boundary condition is implemented based on Eq. (19). In Fig. 5(b), we compare the transient concentration profiles obtained using tDPD with the theoretical solution Eq. (26). An excellent agreement between the tDPD results and the theoretical solution confirms the validity of the boundary method for imposing the correct Neumann boundary condition.

The one-dimensional diffusions validated the tDPD model and boundary methods for implementation of the correct Dirichlet and Neumann boundary conditions. However, these cases are only one-dimensional and do not consider the flow effect, i.e., advection. A more challenging test for the present tDPD model is solving two-dimensional advection-diffusion equations.

Figure 6 illustrates the geometry and boundary conditions of the problem. A fully developed plane Poiseuille flow is considered between two stationary walls. Both walls have no-slip boundary conditions for velocity. For the concentration field, the upper wall has fixed concentration C = 1 while on the lower wall we prescribe the boundary condition dC/dn = 0 for concentration.

Three tDPD simulations with different Peclet numbers Pe = 10, 100, 1000 are performed in a computational domain of 50 × 10 × 10. The tDPD parameters are the same as those used in one-dimensional tests. The kinematic viscosity of the tDPD fluid is 6.62. The wall boundary conditions are imposed based on the boundary method described in Section II B. To generate a fully developed Poiseuille flow, we use periodic boundary condition in x-direction and apply a body force of g_x = 3.50 to drive the fluid, which yields a Reynolds number of Re = U_maxH/ν = 10, where U_max is the centerline velocity and H is the height of the channel. However, for the concentration, we apply non-periodic inflow/outflow boundary conditions in x-direction based on the methodology proposed by Lei et al.³² In practice, to impose a correct inflow boundary condition, we use a buffer region of r_c × 10 × 10 and prescribe C = 0 in the buffer region. The outflow boundary condition is implemented with a maximum buffer length of 2r_c in which dC/dx = 0 is applied. We note that outflow buffer regions with length of 2r_c and r_c yield same results when the Peclet number (Pe = U_maxH/D) is high, e.g., Pe = 100 and 1000, but a shorter buffer region at Pe = 10 results in distortions of the concentration field in the vicinity of the outlet due to stronger diffusion.

The concentration fields obtained by the tDPD simulations at three Peclet numbers are displayed in Figs. 7(a1), 7(b1), and 7(c1). The Reynolds number of the channel flow is fixed at Re = 10, and various Peclet numbers are obtained through the relation Pe = Re ⋅ Sc by varying Schmidt numbers given by Sc = ν/D, in which the diffusion coefficient D is a linear function of Fickian friction coefficient κ shown in Eq. (12).

Advection-diffusion equations are also solved in the same system using spectral element method (SEM),³³ where high-order hierarchical Jacobi polynomials of order 5 have been used for the solution of velocity, pressure, and concentration fields. Comparisons between the results of tDPD and SEM are presented in Figs. 7(a2), 7(b2), and 7(c2). Concentration profiles shown in Figs. 7(a2), 7(b2), and 7(c2) are extracted from the concentration fields along x-direction at three different heights. The results of tDPD simulations indicate very good agreement with SEM results for various Peclet numbers, especially for smaller Peclet numbers. However, at higher Pe = 1000, a very thin boundary layer is formed on the upper wall shown by the concentration field of Fig. 7(c1) resulting in small differences between tDPD and SEM as shown in Fig. 7(c2). The small discrepancy is due to a low number density of particles ρ = 4 in our tDPD simulations that cannot provide enough spatial resolution in the thin layer.

In this section, we present an application of tDPD to a biological system involving many reactant species to demonstrate the capability of the tDPD model. Here, we look at the coagulation cascade, which is one of the major processes in thrombosis. The main product of the cascade is thrombin, which is mainly responsible in transforming fibrinogen to high concentrations of fibrin at the site of injury forming a fibrin network to increase the structural stiffness of thrombus. It is also known that thrombin enhances platelet activations by approximately 20%.³⁴ It should be noted that there are several mathematical models describing the coagulation cascade leading to a system of partial differential equations (PDEs) and ordinary differential equations (ODEs). Among those models, Anand et al.³⁵ proposed a set of 23 coupled ADR equations for describing the evolution of 25 biological reactants involved in a combined model of intrinsic and extrinsic pathways of blood coagulation process and fibrinolysis, which are listed in Table II. We will be focusing on this model since it takes both blood flow and spatial heterogeneity into consideration. We would like to emphasize that, for simplicity, we are looking at blood coagulation cascade in simple Newtonian fluid. Hence, we are not considering the physiological blood flow in small arteries, where non-Newtonian constitutive law for the DPD fluid is needed. Further, a detailed analysis of the model accuracy and effectiveness in capturing the coagulation dynamics is beyond the scope of this study.

The advection-diffusion-reaction equations are in the form of

where C_i denotes the concentration of ith reactant and D_i the corresponding diffusion coefficient. $Q i S$ is the nonlinear source terms representing the production or depletion of C_i due to the enzymatic cascade of reactions. To simplify the formulas in this section, we number the 25 species as shown in Table II and use symbol [i] to denote the concentration of ith reactant.

Equation (27) includes 23 ADR equations involving 25 chemical species. Different species affect each other through chemical reactions. The expressions of the 23 source terms are shown in Table III. The involved reaction kinetics parameters are given in the caption of Table III. The initial conditions of the fibrinogen(I), zymogens (II, V, VIII, IX, X, XI), inhibitors (ATIII, TFPI), and other factors (PC, α₁AT, tPA, PLS, α₂AP) are determined based on normal human physiological values.^15,35,36 We also consider 0.1% initial concentration for the enzymes (IIa, Va, VIIIa, IXa, Xa, XIa), active protein C (APC), and Plasmin (PLA) as shown in Table II in order to initiate the coagulation process. The concentrations of two other chemical species tenase (Z) and prothrombinase (W) are computed through the relations [24] = [3][1]/0.56 and [25] = [5][7]/0.1, respectively.³⁵ 

We consider a flow with maximum velocity of u = 1 mm/s through a micro-channel with width of 20 μm in room temperature. The viscosity of the fluid is 1.0 × 10⁻⁶ m²/s. Taking the diffusivity of the thrombin as 6.47 × 10⁻¹¹ m²/s, the corresponding Peclet and Schmidt numbers are Pe = 309 and Sc = 1.546 × 10⁴. The corresponding Peclet number of various biological species can be computed in a similar way. Further, noting that Pe = Re ⋅ Sc, the value of Reynolds number Re ≈ 0.02 so the flow is approximately in the Stokes regime. The above values are based on the physiologically relevant values of blood velocity and viscosity in arterioles and venules.

As mentioned in Section II A, the effective diffusion coefficient D of a tDPD system is a combined result of the random diffusion D^ξ and the Fickian diffusion D^F. Although the magnitude of D can be linearly adjusted by varying the Fickian friction coefficient κ, D has a minimum of D_min = D^ξ at κ = 0. Therefore, the maximum Schmidt number that a tDPD system can achieve is determined by Sc(κ = 0) = ν/D_min. Table II lists the Schmidt numbers for different chemical species, which are in the order of 10⁴. Such high values of Schmidt numbers in tDPD are achieved through increasing the momentum cutoff radius r_c. We construct a tDPD system in a computational domain of 100 × 10 × 20 with the Schmidt numbers listed in Table II. The parameters of the tDPD interactions are given as ρ = 4.0, k_BT = 1.0, a = 18.75, γ = 4.5, σ = 3.0, r_c = 2.9, r_cc = 1.0, ω_C(r) = (1 − r/r_c), $ω D ( r ) = ω R 2 ( r ) = ( 1 − r / r c ) 0 . 4$⁠, and ω_DC(r) = (1 − r/r_cc)². Note that here we set r_cc < r_c to reduce the computational cost of solving advection-diffusion for 25 species. Computations for kinematic viscosity and coefficient of self diffusion based on this system will result in ν = 138.6 and D = 0.002 61. The maximum Schmidt number is Sc = ν/D = 5.31 × 10⁴, which ensures that this tDPD system can cover all the Schmidt numbers in Table II. The diffusion coefficients for all species can be evaluated by using the method analogous to the reverse Poiseuille flow as described in Section II A, which gives D_i = 0.002 61 + 0.0645κ_i. Since the Schmidt numbers of different species Sc_i are listed in Table II, the diffusion coefficient D_i can be obtained by D_i = ν/Sc_i with the viscosity ν = 138.6. The length and time scales are chosen as L^∗ = 1 μm and t^∗ = 1.386 × 10⁻⁴ s. The typical time scale of blood coagulation dynamics is in the order of minutes, hence making the long-time tDPD computation very expensive and almost impractical. Therefore, in the tDPD modeling, we keep the same Peclet number to capture the correct flow physics, while we accelerate the reaction processes by 1000 times, which compresses the dynamic process of 1 min into 60 ms. This adjustment, however, may reduce the size of the zone in which the reactions occur, and hence the results are qualitative. Although the Peclet number is the relevant nondimensional parameter, we use the physiologic Schmidt and lower Reynolds numbers in the simulations so that the reactions occur at the site of injury before chemical species get advected downstream. Also, the concentrations are scaled by 1.0 nM ⋅ (L^∗)³/ρ with ρ the number density of tDPD particles, and hence a concentration of C^∗ = 1.0 nM is scaled into the concentration unit in our tDPD system. We note that all the parameters without physical units in Tables II–IV are defined in DPD units. The details of mapping physical quantities to tDPD parameters are described in Appendix  B.

Before we perform the simulation of blood coagulation, we estimate the extra computational cost for solving 23 ADR equations by performing a tDPD simulation and a conventional DPD simulation. A small fluid system with size of 30.0 × 33.0 × 10.0 including 39 600 particles is employed to test the computational cost on a single core of Intel Xeon E5-2637 CPU running at 3.50 GHz. We find that the tDPD simulation involving 23 ADR equations (25 chemical species) needs 0.398 s per time step, while the conventional DPD without any chemical species needs 0.192 s per time step. We conclude that the total computational cost of tDPD for solving 23 ADR equations is only doubled compared to the flow solver using the conventional DPD. However, in continuum models, a Poisson solver occupies approximately 60% computation time,²² which makes the additional cost of 23 Poisson solvers over ten times higher than the pure Navier-Stokes solver.

Figure 8 illustrates the initial and boundary conditions for the tDPD simulation of the coagulation and fibrinolysis dynamics in blood flows. The velocity field is periodic in x and z directions with a fully developed Poiseuille flow inside the channel. Concentration of every species at the inlet is set to its initial concentration $C i 0$ and is imposed in a buffer region of length r_c. Further, zero Neumann boundary condition is imposed in a buffer region in the similar way described in Sec. II B. We assume that coagulation occurs at the site of injury (0.1 < x/L_x < 0.2), and the process is initiated and sustained by imposing a boundary flux for 7 species listed in Table IV. Note that the tissue factor (TF-VIIa complex) is one of the important factors triggering the extrinsic pathway of coagulation, whereas other species are involved in the intrinsic (contact) pathway, except for tPA which starts the fibrinolysis (degradation of fibrin) process.

Representative results are shown in Figs. 9 and 10 noting that concentrations are evolving in a 1000 times accelerated time sequence. Figure 9 shows the time evolution contours of fibrin concentration, which can represent the growth of thrombus at the site of injury and its diffusion downstream. Here, we assume a threshold value of 3 μM for fibrin to represent the actual geometry of the clot. More detailed results can be found in Fig. 10 where we plot time-variation of concentrations for thrombin, fibrin, anti-thrombin, and APC at the center of the vessel wall injury. These results are in qualitative agreement with the previous numerical studies.^35,36 The most prominent observation in coagulation dynamics of blood is the significant increase of thrombin concentration in the vicinity of the wall during the initiation phase also known as thrombin burst followed by a drop.³⁶ This burst is also accompanied by a significant increase in fibrin and APC concentrations and a major drop in anti-thrombin concentration. Fibrin concentration shows a slight decrease following the initial peak until it levels off to a steady value, which in turn explains the unchanged concentration contours after Fig. 9(f). All the observations in tDPD simulation with time acceleration are consistent with the results of Bodnár and Sequeira.³⁶