We present a transport dissipative particle dynamics (tDPD) model for simulating mesoscopic problems involving advection-diffusion-reaction (ADR) processes, along with a methodology for implementation of the correct Dirichlet and Neumann boundary conditions in tDPD simulations. tDPD is an extension of the classic dissipative particle dynamics (DPD) framework with extra variables for describing the evolution of concentration fields. The transport of concentration is modeled by a Fickian flux and a random flux between tDPD particles, and the advection is implicitly considered by the movements of these Lagrangian particles. An analytical formula is proposed to relate the tDPD parameters to the effective diffusion coefficient. To validate the present tDPD model and the boundary conditions, we perform three tDPD simulations of one-dimensional diffusion with different boundary conditions, and the results show excellent agreement with the theoretical solutions. We also performed two-dimensional simulations of ADR systems and the tDPD simulations agree well with the results obtained by the spectral element method. Finally, we present an application of the tDPD model to the dynamic process of blood coagulation involving 25 reacting species in order to demonstrate the potential of tDPD in simulating biological dynamics at the mesoscale. We find that the tDPD solution of this comprehensive 25-species coagulation model is only twice as computationally expensive as the conventional DPD simulation of the hydrodynamics only, which is a significant advantage over available continuum solvers.

## I. INTRODUCTION

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.^{1} To this end, mesoscopic methods with stochastic terms are attracting increasing attention as a promising approach for tackling challenging problems in bioengineering and biotechnology.^{2} 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,^{3} and the effects of unsolved degrees of freedom are approximated by stochastic dynamics.^{4} 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.^{5} 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.^{6}

DPD was initially proposed by Hoogerbrugge and Koelman^{3} 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.^{5} By using the Fokker-Planck equation and applying the Mori projection operator, Español^{7} and Marsh *et al.*^{8} 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,^{13} and nanoparticles across lipid membranes.^{14}

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.^{1} 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.^{15} 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.^{15} 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.^{16} 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,^{17} 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.^{4} 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.

. | DPD . | SPH . |
---|---|---|

Major Diff. | Particle stochastic approach | Continuum deterministic approach |

CG force field (bottom-up) governing DPD particles | Discretization of continuum equations (i.e., Navier-Stokes) | |

Inputs | Particle interactions | Equation of state, viscosity, thermal conductivity, etc. |

Outputs | Equation of state, diffusivity, viscosity, thermal | Close to the inputs |

conductivity, etc. | ||

Adv. | No requirements for constitutive equation, which is | Clear physical definition of parameters |

suitable for modeling complex systems and materials | in continuum equations | |

Disadv. | No clear physical definition of parameters; thus, it | Must know the constitutive equation and |

needs mapping DPD units to physical units based on output properties | macroscopic properties given as inputs |

. | DPD . | SPH . |
---|---|---|

Major Diff. | Particle stochastic approach | Continuum deterministic approach |

CG force field (bottom-up) governing DPD particles | Discretization of continuum equations (i.e., Navier-Stokes) | |

Inputs | Particle interactions | Equation of state, viscosity, thermal conductivity, etc. |

Outputs | Equation of state, diffusivity, viscosity, thermal | Close to the inputs |

conductivity, etc. | ||

Adv. | No requirements for constitutive equation, which is | Clear physical definition of parameters |

suitable for modeling complex systems and materials | in continuum equations | |

Disadv. | No clear physical definition of parameters; thus, it | Must know the constitutive equation and |

needs mapping DPD units to physical units based on output properties | macroscopic properties given as inputs |

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.^{22}

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.

## II. GOVERNING EQUATIONS

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^{5}

where *r _{ij}* = |

**r**

_{ij}| with

**r**

_{ij}=

**r**

_{i}−

**r**

_{j},

**e**

_{ij}=

**r**

_{ij}/

*r*is the unit vector from particle

_{ij}*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*is the conservative force parameter,

_{ij}*γ*the dissipative coefficient, and

_{ij}*σ*the strength of random force. Moreover,

_{ij}*C*represents the concentration of one species defined as the number of a chemical species carried by a tDPD particle

_{i}*i*and

*Q*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

_{i}*ρC*where

_{i}*ρ*is the number density of tDPD particles. We note that

*C*can be a vector

_{i}**C**

_{i}containing

*N*components, i.e., {

*C*

_{1},

*C*

_{2}, …,

*C*}

_{N}_{i}when

*N*chemical species are considered.

Based on Fick’s law,^{23} 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

*ϵ*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.

_{ij}^{5}A similar DPD transport model was proposed in Ref. 20 by drawing an analogy with the DPD model of heat conduction.

^{18}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.

^{24}

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*the mass of a single solute molecule. Also,

_{s}*C*and

_{i}*C*are the concentrations on particles

_{j}*i*and

*j*, respectively. Simplified derivations for obtaining Eq. (9) are presented in Appendix A. Since

*m*is small with respect to the mass of a tDPD particle

_{s}*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*becomes comparable to

_{s}*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.^{26} 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.

### A. Mesoscopic friction of concentration

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.*

^{21}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^{5}

where *r _{c}* is the cutoff radius for forces. Also, the macroscopic diffusion coefficient

*D*due to Fickian concentration flux can be computed by

^{F}^{27}

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*less than or equal to

_{cc}*r*. Specifically, when the radial distribution function of ideal gas

_{c}*g*(

*r*) = 1.0 is employed, both

*D*

^{ξ}and

*D*can be evaluated analytically and we have

^{F} where *s*_{1} and *s*_{2} are the exponents of the weighting functions given as *ω _{D}* = (1 −

*r*/

*r*)

_{c}^{s1}and

*ω*= (1 −

_{DC}*r*/

*r*)

_{cc}^{s2}, 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*with temperature and the particle size (i.e.,

^{F}*r*and

_{c}*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*surrounded by other particles with concentration of

_{i}*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*. Then, the variation of concentration on particle

_{cc}*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*∇

^{2}

*C*= −

*Q*, which is the same as the governing equation of Poiseuille flow driven by a body force

^{S}*υ*∇

^{2}

**V**= −

*g*. Therefore, we can use a method analogous to the reverse Poiseuille flow

^{28}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 _{B}T* = 1.0,

*a*= 75

*k*/

_{B}T*ρ*,

*γ*= 4.5,

*σ*

^{2}= 2

*k*,

_{B}Tγ*r*=

_{c}*r*= 1.58,

_{cc}*ω*(

_{C}*r*) = (1 −

*r*/

*r*), $ \omega D ( r ) = \omega R 2 ( r ) = ( 1 \u2212 r / r c ) 0 . 41 $, and

_{c}*ω*(

_{DC}*r*) = (1 −

*r*/

*r*)

_{cc}^{2}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*is applied in the domain of

^{S}*z*> 0 while a concentration sink with same magnitude −

*Q*is applied in the other domain

^{S}*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

^{−3}is the source term,

*C*

_{0}= 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).

### B. Boundary conditions

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.^{29} 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.*^{32} 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.*^{19} 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.

#### 1. Dirichlet boundary condition

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^{19} by

where *C _{i}* is the concentration of particle

*i*,

*C*the expected concentration at the wall surface, and

_{w}*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 _{B}T* = 1.0,

*a*= 18.75,

*γ*= 4.5,

*σ*= 3.0,

*r*= 1.58,

_{c}*ω*(

_{C}*r*) = (1 −

*r*/

*r*), and $ \omega D ( r ) = \omega R 2 ( r ) = ( 1 \u2212 r / r c ) 0 . 41 $, the radial distribution function

_{c}*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}^{2}and

*r*=

_{cc}*r*. We note that though

_{c}*r*is usually set to be

_{cc}*r*for convenience,

_{c}*r*can be different from

_{cc}*r*. The dependence of

_{c}*φ*(

*h*) on the distance

*h*is presented in Fig. 2(b), which is approximated by

*φ*(

*h*) = 4.9562 × 10

^{−2}/

*h*⋅ [1 + 4.1517

*h*/

*r*+ 7.1697(

_{cc}*h*/

*r*)

_{cc}^{2}](1 −

*h*/

*r*)

_{cc}^{4.09}for

*h*<

*r*. As the distance

_{cc}*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.01

*r*) =

_{cc}*φ*(0.01

*r*).

_{cc}#### 2. Neumann boundary condition

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 $ \Phi \xaf ( h ) =\Phi ( h ) / \u222b 0 r c c \Phi ( x ) dx$. Then, the integral of $ \Phi \xaf ( h ) $ is equal to one. Using the computed *g*(*r*) shown in Fig. 2(a) and *ω _{DC}*(

*r*) = (1 −

*r*/

*r*)

_{cc}^{2}, the function $ \Phi \xaf ( h ) $ can be obtained through Eq. (18), as shown in Fig. 2(b). The computed $ \Phi \xaf ( h ) $ is approximated by $ \Phi \xaf ( h ) = \u2211 i = 0 6 p i ( 1 \u2212 h / r c c ) 7 \u2212 i $ with the coefficients

**p**= {22.745 45, − 62.6256, 55.1365, − 19.2584, 6.2662, − 0.502 42, 0.0101} for

*h*≤

*r*and $ \Phi \xaf ( h ) =0$ for

_{cc}*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*is distributed onto the fluid particles in the vicinity of the wall weighted by $ \Phi \xaf ( h ) $. Let

^{W}*ρ*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 ) = \u2211 i = 0 6 p i ( 1 \u2212 h / r c ) 7 \u2212 i $ with the coefficients **p** = {525.6666, − 1457.0412, 1485.7229, − 927.0881, 453.3782, − 7.2616, 0.177 58},

where *v _{n}*

**e**

_{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*

_{n}**e**

_{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.6212

*h*/

*r*+ 2.2090(

_{c}*h*/

*r*)

_{c}^{2}](1 −

*h*/

*r*)

_{c}^{3.45}and

*γ*(

_{n}*h*) = 20.7789/

*h*⋅ [1 + 2.1245

*h*/

*r*+ 5.4566(

_{c}*h*/

*r*)

_{c}^{2}](1 −

*h*/

*r*)

_{c}^{2.30}for

*h*<

*r*, respectively.

_{c}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 $ \sigma \tau 2 ( h ) =2 k B T \gamma \tau ( h ) $ and $ \sigma n 2 ( h ) =2 k B T \gamma 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.05

*r*) and

_{c}*γ*(0.1

_{n}*r*) to avoid instability of tDPD simulations induced by some extremely large random numbers.

_{c}The boundary method based on the effective boundary forces has been well validated by Lei *et al.*^{32} and Li *et al.*^{19} 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.

## III. RESULTS

### A. Validation for one-dimensional diffusion

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*_{1} and *C*_{2} are the concentration of species 1 and 2, respectively. *D*_{11} and *D*_{22} 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*_{11} and *D*_{22} 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} = 1.0 and *C*_{2} = 0.0 in the domain of *x* < 0 and *C*_{1} = 0.0 and *C*_{2} = 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 _{B}T* = 1.0,

*a*= 75

*k*/

_{B}T*ρ*,

*γ*= 4.5,

*σ*

^{2}= 2

*k*,

_{B}Tγ*r*=

_{c}*r*= 1.58,

_{cc}*ω*(

_{C}*r*) = (1 −

*r*/

*r*), $ \omega D ( r ) = \omega R 2 ( r ) = ( 1 \u2212 r / r c ) 0 . 41 $, and

_{c}*ω*(

_{DC}*r*) = (1 −

*r*/

*r*)

_{cc}^{2}. Unless otherwise specified, all the simulations in the present paper are based on this set of parameters. Also,

*κ*

_{11}= 0.04 results in

*D*

_{11}= 0.05 and

*κ*

_{22}= 0.62 yields

*D*

_{22}= 0.50. The velocity-Verlet algorithm

^{5}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/ \pi \u22c5 \u222b 0 x e \u2212 \tau 2 d\tau $ 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*_{0}, an analytical solution for the transient concentration profile can be obtained,

where *β _{n}* =

*nπ*/

*L*with

_{x}*L*being the length of the computational domain in the

_{x}*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π*/(2

*L*) with

_{x}*L*being the length of the computational domain in the

_{x}*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.

### B. Validation for two-dimensional advection-diffusion

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*/

_{max}H*ν*= 10, where

*U*is the centerline velocity and

_{max}*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.*

^{32}In practice, to impose a correct inflow boundary condition, we use a buffer region of

*r*× 10 × 10 and prescribe

_{c}*C*= 0 in the buffer region. The outflow boundary condition is implemented with a maximum buffer length of 2

*r*in which

_{c}*dC*/

*dx*= 0 is applied. We note that outflow buffer regions with length of 2

*r*and

_{c}*r*yield same results when the Peclet number (

_{c}*Pe*=

*U*/

_{max}H*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),^{33} 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.

### C. Application to biochemical systems

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%.^{34} 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.*^{35} 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.

Index . | Species . | Name . | $ C i 0 $ . | Sc
. |
---|---|---|---|---|

1 | IXa | Enzym IXa | 0.09 | 15 949 |

2 | IX | Zymogen IX | 90 | 17 762 |

3 | VIIIa | Enzym VIIIa | 0.0007 | 25 510 |

4 | VIII | Zymogen VIII | 0.7 | 32 051 |

5 | Va | Enzym Va | 0.02 | 26 178 |

6 | V | Zymogen V | 20 | 17 762 |

7 | Xa | Enzym Xa | 0.17 | 13 568 |

8 | X | Zymogen X | 170 | 17 762 |

9 | IIa | Thrombin | 1.4 | 15 456 |

10 | II | Prothrombin | 1 400 | 19 194 |

11 | Ia | Fibrin | 7.0 | 40 486 |

12 | I | Fibrinogen | 7 000 | 32 258 |

13 | XIa | Enzym XIa | 0.03 | 20 000 |

14 | XI | Zymogen XI | 30 | 25 189 |

15 | ATIII | AntiThrombin-III | 2 410 | 17 953 |

16 | TFPI | Tissue factor pathway inhibitor | 2.5 | 15 873 |

17 | APC | Active protein C | 0.06 | 18 182 |

18 | PC | Protein C | 60 | 18 382 |

19 | α_{1} AT | α_{1}-AntiTrypsin | 45 000 | 17 182 |

20 | tPA | Tissue plasminogen activator | 0.08 | 18 939 |

21 | PLA | Plasmin | 2.18 | 20 284 |

22 | PLS | Plasminogen | 2 180 | 20 790 |

23 | α_{2} AP | α_{2}-AntiPlasmin | 105 | 19 048 |

24 | Z | Tenase | 1.125 × 10^{−4} | … |

25 | W | Prothrombinase | 0.034 | … |

Index . | Species . | Name . | $ C i 0 $ . | Sc
. |
---|---|---|---|---|

1 | IXa | Enzym IXa | 0.09 | 15 949 |

2 | IX | Zymogen IX | 90 | 17 762 |

3 | VIIIa | Enzym VIIIa | 0.0007 | 25 510 |

4 | VIII | Zymogen VIII | 0.7 | 32 051 |

5 | Va | Enzym Va | 0.02 | 26 178 |

6 | V | Zymogen V | 20 | 17 762 |

7 | Xa | Enzym Xa | 0.17 | 13 568 |

8 | X | Zymogen X | 170 | 17 762 |

9 | IIa | Thrombin | 1.4 | 15 456 |

10 | II | Prothrombin | 1 400 | 19 194 |

11 | Ia | Fibrin | 7.0 | 40 486 |

12 | I | Fibrinogen | 7 000 | 32 258 |

13 | XIa | Enzym XIa | 0.03 | 20 000 |

14 | XI | Zymogen XI | 30 | 25 189 |

15 | ATIII | AntiThrombin-III | 2 410 | 17 953 |

16 | TFPI | Tissue factor pathway inhibitor | 2.5 | 15 873 |

17 | APC | Active protein C | 0.06 | 18 182 |

18 | PC | Protein C | 60 | 18 382 |

19 | α_{1} AT | α_{1}-AntiTrypsin | 45 000 | 17 182 |

20 | tPA | Tissue plasminogen activator | 0.08 | 18 939 |

21 | PLA | Plasmin | 2.18 | 20 284 |

22 | PLS | Plasminogen | 2 180 | 20 790 |

23 | α_{2} AP | α_{2}-AntiPlasmin | 105 | 19 048 |

24 | Z | Tenase | 1.125 × 10^{−4} | … |

25 | W | Prothrombinase | 0.034 | … |

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

where *C _{i}* denotes the concentration of

*i*th reactant and

*D*the corresponding diffusion coefficient. $ Q i S $ is the nonlinear source terms representing the production or depletion of

_{i}*C*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}*i*] to denote the concentration of

*i*th 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, *α*_{1}AT, tPA, PLS, *α*_{2}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.^{35}

$ Q 1 S = ( k 9 [ 13 ] [ 2 ] ) / ( K 9 M + [ 2 ] ) \u2212 h 9 [ 1 ] [ 15 ] $ |

$ Q 2 S =\u2212 ( k 9 [ 13 ] [ 2 ] ) / ( K 9 M + [ 2 ] ) $ |

$ Q 3 S = ( k 8 [ 9 ] [ 4 ] ) / ( K 8 M + [ 4 ] ) \u2212 h 8 [ 3 ] \u2212 ( h C 8 [ 17 ] [ 3 ] ) / ( H C 8 M + [ 3 ] ) $ |

$ Q 4 S = ( k 8 [ 9 ] [ 4 ] ) / ( K 8 M + [ 4 ] ) $ |

$ Q 5 S = ( k 5 [ 9 ] [ 6 ] ) / ( K 5 M + [ 6 ] ) \u2212 h 5 [ 5 ] \u2212 ( h C 5 [ 17 ] [ 5 ] ) / ( H C 5 M + [ 5 ] ) $ |

$ Q 6 S =\u2212 ( k 5 [ 9 ] [ 6 ] ) / ( K 5 M + [ 6 ] ) $ |

$ Q 7 S = ( k 10 [ 24 ] [ 8 ] ) / ( K 10 M + [ 8 ] ) \u2212 h 10 [ 7 ] [ 15 ] \u2212 h TFPI [ 16 ] [ 7 ] $ |

$ Q 8 S =\u2212 ( k 10 [ 24 ] [ 8 ] ) / ( K 10 M + [ 8 ] ) $ |

$ Q 9 S = ( k 2 [ 25 ] [ 10 ] ) / ( K 2 M + [ 10 ] ) \u2212 h 2 [ 9 ] [ 15 ] $ |

$ Q 10 S =\u2212 ( k 2 [ 25 ] [ 10 ] ) / ( K 2 M + [ 10 ] ) $ |

$ Q 11 S = ( k 1 [ 9 ] [ 12 ] ) / ( K 1 M + [ 12 ] ) \u2212 ( h 1 [ 21 ] [ 11 ] ) / ( H 1 M + [ 11 ] ) $ |

$ Q 12 S =\u2212 ( k 1 [ 9 ] [ 12 ] ) / ( K 1 M + [ 12 ] ) $ |

$ Q 13 S = ( k 11 [ 9 ] [ 14 ] ) / ( K 11 M + [ 14 ] ) \u2212 h 11 A 3 [ 13 ] [ 15 ] \u2212 h 11 L 1 [ 13 ] [ 19 ] $ |

$ Q 14 S =\u2212 ( k 11 [ 9 ] [ 14 ] ) / ( K 11 M + [ 14 ] ) $ |

$ Q 15 S =\u2212 ( h 9 [ 1 ] + h 10 [ 7 ] + h 2 [ 9 ] + h 11 A 3 [ 13 ] ) [ 15 ] $ |

$ Q 16 S =\u2212 h TFPI [ 16 ] [ 7 ] $ |

$ Q 17 S = ( k PC [ 9 ] [ 18 ] ) / ( K PCM + [ 18 ] ) \u2212 h PC [ 17 ] [ 19 ] $ |

$ Q 18 S =\u2212 ( k PC [ 9 ] [ 18 ] ) / ( K PCM + [ 18 ] ) $ |

$ Q 19 S =\u2212 h PC [ 17 ] [ 19 ] \u2212 h 11 L 1 [ 13 ] [ 19 ] $ |

$ Q 20 S =0$ |

$ Q 21 S = ( k PLA [ 20 ] [ 22 ] ) / ( K PLAM + [ 22 ] ) \u2212 h PLA [ 21 ] [ 23 ] $ |

$ Q 22 S =\u2212 ( k PLA [ 20 ] [ 22 ] ) / ( K PLAM + [ 22 ] ) $ |

$ Q 23 S =\u2212 h PLA [ 21 ] [ 23 ] $ |

$ Q 1 S = ( k 9 [ 13 ] [ 2 ] ) / ( K 9 M + [ 2 ] ) \u2212 h 9 [ 1 ] [ 15 ] $ |

$ Q 2 S =\u2212 ( k 9 [ 13 ] [ 2 ] ) / ( K 9 M + [ 2 ] ) $ |

$ Q 3 S = ( k 8 [ 9 ] [ 4 ] ) / ( K 8 M + [ 4 ] ) \u2212 h 8 [ 3 ] \u2212 ( h C 8 [ 17 ] [ 3 ] ) / ( H C 8 M + [ 3 ] ) $ |

$ Q 4 S = ( k 8 [ 9 ] [ 4 ] ) / ( K 8 M + [ 4 ] ) $ |

$ Q 5 S = ( k 5 [ 9 ] [ 6 ] ) / ( K 5 M + [ 6 ] ) \u2212 h 5 [ 5 ] \u2212 ( h C 5 [ 17 ] [ 5 ] ) / ( H C 5 M + [ 5 ] ) $ |

$ Q 6 S =\u2212 ( k 5 [ 9 ] [ 6 ] ) / ( K 5 M + [ 6 ] ) $ |

$ Q 7 S = ( k 10 [ 24 ] [ 8 ] ) / ( K 10 M + [ 8 ] ) \u2212 h 10 [ 7 ] [ 15 ] \u2212 h TFPI [ 16 ] [ 7 ] $ |

$ Q 8 S =\u2212 ( k 10 [ 24 ] [ 8 ] ) / ( K 10 M + [ 8 ] ) $ |

$ Q 9 S = ( k 2 [ 25 ] [ 10 ] ) / ( K 2 M + [ 10 ] ) \u2212 h 2 [ 9 ] [ 15 ] $ |

$ Q 10 S =\u2212 ( k 2 [ 25 ] [ 10 ] ) / ( K 2 M + [ 10 ] ) $ |

$ Q 11 S = ( k 1 [ 9 ] [ 12 ] ) / ( K 1 M + [ 12 ] ) \u2212 ( h 1 [ 21 ] [ 11 ] ) / ( H 1 M + [ 11 ] ) $ |

$ Q 12 S =\u2212 ( k 1 [ 9 ] [ 12 ] ) / ( K 1 M + [ 12 ] ) $ |

$ Q 13 S = ( k 11 [ 9 ] [ 14 ] ) / ( K 11 M + [ 14 ] ) \u2212 h 11 A 3 [ 13 ] [ 15 ] \u2212 h 11 L 1 [ 13 ] [ 19 ] $ |

$ Q 14 S =\u2212 ( k 11 [ 9 ] [ 14 ] ) / ( K 11 M + [ 14 ] ) $ |

$ Q 15 S =\u2212 ( h 9 [ 1 ] + h 10 [ 7 ] + h 2 [ 9 ] + h 11 A 3 [ 13 ] ) [ 15 ] $ |

$ Q 16 S =\u2212 h TFPI [ 16 ] [ 7 ] $ |

$ Q 17 S = ( k PC [ 9 ] [ 18 ] ) / ( K PCM + [ 18 ] ) \u2212 h PC [ 17 ] [ 19 ] $ |

$ Q 18 S =\u2212 ( k PC [ 9 ] [ 18 ] ) / ( K PCM + [ 18 ] ) $ |

$ Q 19 S =\u2212 h PC [ 17 ] [ 19 ] \u2212 h 11 L 1 [ 13 ] [ 19 ] $ |

$ Q 20 S =0$ |

$ Q 21 S = ( k PLA [ 20 ] [ 22 ] ) / ( K PLAM + [ 22 ] ) \u2212 h PLA [ 21 ] [ 23 ] $ |

$ Q 22 S =\u2212 ( k PLA [ 20 ] [ 22 ] ) / ( K PLAM + [ 22 ] ) $ |

$ Q 23 S =\u2212 h PLA [ 21 ] [ 23 ] $ |

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^{−6} m^{2}/s. Taking the diffusivity of the thrombin as 6.47 × 10^{−11} m^{2}/s, the corresponding Peclet and Schmidt numbers are *Pe* = 309 and *Sc* = 1.546 × 10^{4}. 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*. Table II lists the Schmidt numbers for different chemical species, which are in the order of 10

_{min}^{4}. Such high values of Schmidt numbers in tDPD are achieved through increasing the momentum cutoff radius

*r*. 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

_{c}*ρ*= 4.0,

*k*= 1.0,

_{B}T*a*= 18.75,

*γ*= 4.5,

*σ*= 3.0,

*r*= 2.9,

_{c}*r*= 1.0,

_{cc}*ω*(

_{C}*r*) = (1 −

*r*/

*r*), $ \omega D ( r ) = \omega R 2 ( r ) = ( 1 \u2212 r / r c ) 0 . 4 $, and

_{c}*ω*(

_{DC}*r*) = (1 −

*r*/

*r*)

_{cc}^{2}. Note that here we set

*r*<

_{cc}*r*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

_{c}*ν*= 138.6 and

*D*= 0.002 61. The maximum Schmidt number is

*Sc*=

*ν*/

*D*= 5.31 × 10

^{4}, 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*= 0.002 61 + 0.0645

_{i}*κ*. Since the Schmidt numbers of different species

_{i}*Sc*are listed in Table II, the diffusion coefficient

_{i}*D*can be obtained by

_{i}*D*=

_{i}*ν*/

*Sc*with the viscosity

_{i}*ν*= 138.6. The length and time scales are chosen as

*L*

^{∗}= 1

*μ*m and

*t*

^{∗}= 1.386 × 10

^{−4}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*

^{∗})

^{3}/

*ρ*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.

Index j
. | Boundary flux terms Q
. _{j} |
---|---|

1 | (k_{7,9}[2]^{W}[TF-VIIa]^{W})/(K_{7,9M} + [2]^{W}) L |

2 | −(k_{7,9}[2]^{W}[TF-VIIa]^{W})/(K_{7,9M} + [2]^{W}) L |

7 | (k_{7,10}[8]^{W}[TF-VIIa]^{W})/(K_{7,10M} + [8]^{W}) L |

8 | −(k_{7,10}[8]^{W}[TF-VIIa]^{W})/(K_{7,10M} + [8]^{W}) L |

13 | (ϕ_{11}[14]^{W}[XIIa]^{W})/(Φ_{11M} + [14]^{W}) L |

14 | −(ϕ_{11}[14]^{W}[XIIa]^{W})/(Φ_{11M} + [14]^{W}) L |

20 | $ ( k tPA C + k tPA IIa [ 9 ] W + k tPA Ia [ 11 ] W ) [ ENDO ] W L$ |

Index j
. | Boundary flux terms Q
. _{j} |
---|---|

1 | (k_{7,9}[2]^{W}[TF-VIIa]^{W})/(K_{7,9M} + [2]^{W}) L |

2 | −(k_{7,9}[2]^{W}[TF-VIIa]^{W})/(K_{7,9M} + [2]^{W}) L |

7 | (k_{7,10}[8]^{W}[TF-VIIa]^{W})/(K_{7,10M} + [8]^{W}) L |

8 | −(k_{7,10}[8]^{W}[TF-VIIa]^{W})/(K_{7,10M} + [8]^{W}) L |

13 | (ϕ_{11}[14]^{W}[XIIa]^{W})/(Φ_{11M} + [14]^{W}) L |

14 | −(ϕ_{11}[14]^{W}[XIIa]^{W})/(Φ_{11M} + [14]^{W}) L |

20 | $ ( k tPA C + k tPA IIa [ 9 ] W + k tPA Ia [ 11 ] W ) [ ENDO ] W L$ |

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,^{22} 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*< 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.

_{x}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.^{36} 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.^{36}

It is worth noting that the treatment of “time acceleration” will result in quantitative error. However, the main purpose of this case is to demonstrate how to apply the tDPD model to complex diffusion-reaction systems. With “time acceleration”, the designed objective of Section III C has been achieved. The error induced by the “time acceleration” can be eliminated by realtime simulation using large-scale parallel computations, which are beyond the scope of the present work.

## IV. SUMMARY AND DISCUSSION

A tDPD model for simulating mesoscopic fluid systems involving ADR processes has been developed. Each tDPD particle is associated with extra variables for carrying concentrations to describe the evolution of concentration fields. The exchange of concentration between tDPD particles is modeled by a Fickian flux and a random flux, and an analytical formula is provided to relate the mesoscopic parameters to the effective diffusion coefficient *D*, which is tunable by linearly varying the Fickian friction coefficient *κ*.

A methodology for imposing the correct Dirichlet and Neumann boundary conditions has been developed. In particular, the effects of a real solid wall are replaced by effective boundary forces and concentration fluxes, which are extracted from the fluid-solid particle interactions using the continuum approximation. Simulations of one-dimensional diffusion with different boundary conditions were carried out to validate the present tDPD model and the boundary methods for imposing the Dirichlet and Neumann boundary conditions. Furthermore, the tDPD model was used to solve two-dimensional advection-diffusion equations with various Peclet numbers, where tDPD results showed very good agreement with the results obtained by the spectral element method.

Since the particle-based tDPD method satisfies the conservation of concentration automatically, it provides an economical way to solve ADR equations with a large number of species. An application of the tDPD model to the dynamic process of blood coagulation was performed to demonstrate the promising biological applications of tDPD. Results revealed that the tDPD simulation provided qualitatively correct evolution of fibrin concentration initialized by a injured vessel wall in flowing blood. Moreover, the tDPD simulation confirmed that the total computational cost for solving 23 ADR equations is only doubled compared to the flow solver, unlike the continuum model requiring 23 Poisson/Helmholtz solvers making the computational cost over ten times higher than the Navier-Stokes solver.

Although we demonstrated the implementation of tDPD simulation with uncoupled linear diffusion equations in which the diffusion coefficients *D* are constant and the Fickian friction coefficient *κ* is a diagonal matrix, tDPD can be extended to modeling nonlinear systems, where diffusivity is a function of concentration, position, and time, i.e., *D* = *D*(*C*, **x**, *t*), because of the linear dependence of *D* on the tDPD parameter *κ*. Furthermore, as future works, it is also interesting to consider a dense mixture involving multiple cross terms and chemical interactions between different species, where the Fickian friction coefficient *κ* becomes an *N* × *N* matrix to represent interdiffusivities, and there is coupling between various chemical species.

## Acknowledgments

This work was primarily supported by NIH (Grant No. 1U01HL116323-01) and the DOE Center on Mathematics for Mesoscopic Modeling of Materials (CM4). Computational resources were provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program and TACC/STAMPEDE through the XSEDE Grant (Grant No. TG-DMS140007). Z. Li would like to acknowledge helpful discussions with Dr. Wenxiao Pan and Professor Bruce Caswell. A. Yazdani would like to thank Dr. Hessam Babaee for the support he gave for the spectral element solver.

### APPENDIX A: CORRELATION OF CONCENTRATION FLUCTUATIONS

In this appendix, we present simplified derivations for obtaining the concentration fluctuations given by Eqs. (7) and (9). According to fluctuating hydrodynamics, the fluctuations can be described by random noise terms whose correlation functions are determined by a fluctuation-dissipation theorem.^{26}

By applying the local-equilibrium assumption to the classical thermodynamics Gibbs-Duhem equation for the specific entropy *s*, we have the balance of entropy given by^{26}

where *T* is the temperature, *u* the energy, *p* the pressure, and $v= \rho 0 \u2212 1 $ the specific volume. Also, *μ _{k}* are the chemical potentials defined per unit mass, and

*ω*the corresponding concentrations in terms of mass fraction. Considering the balance laws in fluids together with phenomenological relations, the correlation of random concentration fluxes

_{k}*Q*is given by the fluctuation-dissipation theorem in the form of

^{R} where *D* is the diffusion coefficient. For more details on the derivations leading to Eq. (A2) from Eq. (A1), we refer interested readers to a book by Ortiz de Zárate and Sengers.^{26}

Considering dilute solutions, the chemical potential of a species can be represented by

where *x* is the molecular fraction. We note that $ m s \u2212 1 $ appears in Eq. (A3) because the chemical potential has been defined per unit mass, and *m _{s}* is the mass of a solute molecule. Let

*c*be the molar concentration in terms of molecules (or moles) per unit volume, then

*c*is related to the mass fraction

*ω*by

where *ρ*_{0} is the density of the solution. By assuming that the mole fraction of the solvent in dilute solutions can be approximated by unity, then we have *x* ≈ *c*. Substituting Eqs. (A3) and (A4) into Eq. (A2) results in

We note that the correlation of concentration fluctuation given by Eq. (A5) includes the equilibrium concentration *c*. Therefore, it is only valid for fluid in thermodynamic equilibrium or with the approximation of locally thermodynamic equilibrium. Moreover, the fluctuation-dissipation theorem is valid only for large cell size limit in which fluctuations are small. To this end, Eq. (A5) cannot be used for giant stochastic fluxes.

In the tDPD model, the diffusion coefficient *D* in Eq. (A5) is replaced by the strength of Fickian friction coefficient and the volume concentration becomes *c* = *ρC* with *ρ* the number density of tDPD particles. Therefore, we have the fluctuation-dissipation relationship for concentration fluxes in tDPD given by

### APPENDIX B: MAPPING OF UNITS

Appendix B presents the details of mapping physical quantities to tDPD parameters that are used in Section III C. Table V lists both the physical units and the tDPD units of the parameters, which appear in Tables III and IV.

. | Physical units . | tDPD units (accelerated) . |
---|---|---|

k_{9} | 11 min^{−1} | 2.54 × 10^{−5} (2.54 × 10^{−2}) |

K_{9M} | 160 nM | 160 |

h_{9} | 0.0162 nM^{−1} min^{−1} | 3.74 × 10^{−8} (3.74 × 10^{−5}) |

k_{8} | 194.4 min^{−1} | 4.49 × 10^{−4} (0.449) |

K_{8M} | 1.12 × 10^{5} nM | 1.12 × 10^{5} |

h_{8} | 0.222 min^{−1} | 5.13 × 10^{−7} (5.13 × 10^{−4}) |

h_{C8} | 10.2 min^{−1} | 2.36 × 10^{−5} (2.36 × 10^{−2}) |

H_{C8M} | 14.6 nM | 14.6 |

k_{5} | 27.0 min^{−1} | 6.24 × 10^{−5} (6.24 × 10^{−2}) |

K_{5M} | 140.5 nM | 140.5 |

h_{5} | 0.17 min^{−1} | 3.93 × 10^{−7} (3.93 × 10^{−4}) |

h_{C5} | 10.2 min^{−1} | 2.36 × 10^{−5} (2.36 × 10^{−2}) |

H_{C5M} | 14.6 nM | 14.6 |

k_{10} | 2391 min^{−1} | 5.523 × 10^{−3} (5.523) |

K_{10M} | 160 nM | 160 |

h_{10} | 0.347 nM^{−1} min^{−1} | 8.01 × 10^{−7} (8.01 × 10^{−4}) |

h _{TFPI} | 0.48 nM^{−1} min^{−1} | 1.11 × 10^{−6} (1.11 × 10^{−3}) |

k_{2} | 1344 min^{−1} | 3.105 × 10^{−3} (3.105) |

K_{2M} | 1060 nM | 1060 |

h_{2} | 0.714 nM^{−1} min^{−1} | 1.65 × 10^{−6} (1.65 × 10^{−3}) |

k_{1} | 3540 min^{−1} | 8.177 × 10^{−3} (8.177) |

K_{1M} | 3160 nM | 3160 |

h_{1} | 1500 min^{−1} | 3.456 × 10^{−3} (3.456) |

K_{1M} | 2.50 × 10^{5} nM | 2.50 × 10^{5} |

k_{11} | 0.0078 min^{−1} | 1.80 × 10^{−8} (1.80 × 10^{−5}) |

K_{11M} | 50 nM | 50 |

$ h 11 A 3 $ | 1.6 × 10^{−3} nM^{−1} min^{−1} | 3.70 × 10^{−9} (3.70 × 10^{−6}) |

$ h 11 L 1 $ | 1.3 × 10^{−5} nM^{−1} min^{−1} | 3.00 × 10^{−11} (3.00 × 10^{−8}) |

k _{PC} | 39.0 min^{−1} | 9.01 × 10^{−5} (9.01 × 10^{−2}) |

K _{PCM} | 3190 nM | 3190 |

h _{PC} | 6.6 × 10^{−7} nM^{−1} min^{−1} | 1.52 × 10^{−12} (1.52 × 10^{−9}) |

k _{PLA} | 12.0 min^{−1} | 2.77 × 10^{−5} (2.77 × 10^{−2}) |

K _{PLAM} | 18.0 nM | 18.0 |

h _{PLA} | 0.096 nM^{−1} min^{−1} | 2.22 × 10^{−7} (2.22 × 10^{−4}) |

k_{7,9} | 32.4 min^{−1} | 7.48 × 10^{−5} (7.48 × 10^{−2}) |

K_{7,9M} | 24.0 nM | 24.0 |

k_{7,10} | 103 min^{−1} | 2.38 × 10^{−4} (0.238) |

K_{7,10M} | 240.0 nM | 240.0 |

ϕ_{11} | 0.034 min^{−1} | 7.85 × 10^{−8} (7.85 × 10^{−5}) |

Φ_{11M} | 2000 nM | 2000 |

$ k tPA C $ | 6.52 × 10^{−13} nM m^{2} min^{−1} | 1.51 × 10^{−6} (1.51 × 10^{−3}) |

$ k tPA IIa $ | 9.27 × 10^{−12}e^{−134.8(t−t0)} | 2.14 × 10^{−5}e^{−972 580(t−t0)} |

m^{2} min^{−1} | (2.14 × 10^{−2}e^{−972.58(t−t0)}) | |

$ k tPA Ia $ | 5.059 × 10^{−18} m^{2} min^{−1} | 1.17 × 10^{−11} (1.17 × 10^{−8}) |

[TF-VIIa] | 0.25 nM | 0.25 |

[ENDO] | 2.0 × 10^{9} cells/m^{2} | 2.0 × 10^{−3} |

[XIIa] | 375 nM | 375 |

. | Physical units . | tDPD units (accelerated) . |
---|---|---|

k_{9} | 11 min^{−1} | 2.54 × 10^{−5} (2.54 × 10^{−2}) |

K_{9M} | 160 nM | 160 |

h_{9} | 0.0162 nM^{−1} min^{−1} | 3.74 × 10^{−8} (3.74 × 10^{−5}) |

k_{8} | 194.4 min^{−1} | 4.49 × 10^{−4} (0.449) |

K_{8M} | 1.12 × 10^{5} nM | 1.12 × 10^{5} |

h_{8} | 0.222 min^{−1} | 5.13 × 10^{−7} (5.13 × 10^{−4}) |

h_{C8} | 10.2 min^{−1} | 2.36 × 10^{−5} (2.36 × 10^{−2}) |

H_{C8M} | 14.6 nM | 14.6 |

k_{5} | 27.0 min^{−1} | 6.24 × 10^{−5} (6.24 × 10^{−2}) |

K_{5M} | 140.5 nM | 140.5 |

h_{5} | 0.17 min^{−1} | 3.93 × 10^{−7} (3.93 × 10^{−4}) |

h_{C5} | 10.2 min^{−1} | 2.36 × 10^{−5} (2.36 × 10^{−2}) |

H_{C5M} | 14.6 nM | 14.6 |

k_{10} | 2391 min^{−1} | 5.523 × 10^{−3} (5.523) |

K_{10M} | 160 nM | 160 |

h_{10} | 0.347 nM^{−1} min^{−1} | 8.01 × 10^{−7} (8.01 × 10^{−4}) |

h _{TFPI} | 0.48 nM^{−1} min^{−1} | 1.11 × 10^{−6} (1.11 × 10^{−3}) |

k_{2} | 1344 min^{−1} | 3.105 × 10^{−3} (3.105) |

K_{2M} | 1060 nM | 1060 |

h_{2} | 0.714 nM^{−1} min^{−1} | 1.65 × 10^{−6} (1.65 × 10^{−3}) |

k_{1} | 3540 min^{−1} | 8.177 × 10^{−3} (8.177) |

K_{1M} | 3160 nM | 3160 |

h_{1} | 1500 min^{−1} | 3.456 × 10^{−3} (3.456) |

K_{1M} | 2.50 × 10^{5} nM | 2.50 × 10^{5} |

k_{11} | 0.0078 min^{−1} | 1.80 × 10^{−8} (1.80 × 10^{−5}) |

K_{11M} | 50 nM | 50 |

$ h 11 A 3 $ | 1.6 × 10^{−3} nM^{−1} min^{−1} | 3.70 × 10^{−9} (3.70 × 10^{−6}) |

$ h 11 L 1 $ | 1.3 × 10^{−5} nM^{−1} min^{−1} | 3.00 × 10^{−11} (3.00 × 10^{−8}) |

k _{PC} | 39.0 min^{−1} | 9.01 × 10^{−5} (9.01 × 10^{−2}) |

K _{PCM} | 3190 nM | 3190 |

h _{PC} | 6.6 × 10^{−7} nM^{−1} min^{−1} | 1.52 × 10^{−12} (1.52 × 10^{−9}) |

k _{PLA} | 12.0 min^{−1} | 2.77 × 10^{−5} (2.77 × 10^{−2}) |

K _{PLAM} | 18.0 nM | 18.0 |

h _{PLA} | 0.096 nM^{−1} min^{−1} | 2.22 × 10^{−7} (2.22 × 10^{−4}) |

k_{7,9} | 32.4 min^{−1} | 7.48 × 10^{−5} (7.48 × 10^{−2}) |

K_{7,9M} | 24.0 nM | 24.0 |

k_{7,10} | 103 min^{−1} | 2.38 × 10^{−4} (0.238) |

K_{7,10M} | 240.0 nM | 240.0 |

ϕ_{11} | 0.034 min^{−1} | 7.85 × 10^{−8} (7.85 × 10^{−5}) |

Φ_{11M} | 2000 nM | 2000 |

$ k tPA C $ | 6.52 × 10^{−13} nM m^{2} min^{−1} | 1.51 × 10^{−6} (1.51 × 10^{−3}) |

$ k tPA IIa $ | 9.27 × 10^{−12}e^{−134.8(t−t0)} | 2.14 × 10^{−5}e^{−972 580(t−t0)} |

m^{2} min^{−1} | (2.14 × 10^{−2}e^{−972.58(t−t0)}) | |

$ k tPA Ia $ | 5.059 × 10^{−18} m^{2} min^{−1} | 1.17 × 10^{−11} (1.17 × 10^{−8}) |

[TF-VIIa] | 0.25 nM | 0.25 |

[ENDO] | 2.0 × 10^{9} cells/m^{2} | 2.0 × 10^{−3} |

[XIIa] | 375 nM | 375 |

The length scale is determined by the size of physical system. In Section III C, a channel with length *L _{x}* = 100

*μ*m and height

*L*= 20

_{z}*μ*m is considered and the length unit is chosen as

*L*

^{∗}= 1

*μ*m. Therefore, a tDPD system is constructed in a computational domain of 100 × 10 × 20 to represent the physical system.

The time unit of tDPD system is determined by

where *ν*_{0} = 1.0 × 10^{−6} m^{2}/s is the kinematic viscosity of fluid in physical units. Since the length unit *L*^{∗} = 1 *μ*m and the kinematic viscosity of the tDPD system is *ν _{DPD}* = 138.6, we have the time unit

*t*

^{∗}= 1.386 × 10

^{−4}s.

Then, all the physical units of length can be scaled by *L*^{∗} and the physical units of time can be scaled by *t*^{∗} to obtain corresponding tDPD units. For example, 1 m/*L*^{∗} = 1 m/1.0 × 10^{−6} m = 1.0 × 10^{6}, 1 m^{3}/(*L*^{∗})^{3} = 1.0 × 10^{18}, 1 min/*t*^{∗} = 60 s/1.386 × 10^{−4} s = 4.33 × 10^{5}, and 1 min^{−1} ⋅ *t*^{∗} = 1.386 × 10^{−4} s/60 s = 2.31 × 10^{−6}.

Moreover, the volume concentration unit is chosen as *C*^{∗} = 1 nM, which indicates that the number of molecules per tDPD particle should be scaled by *n*^{∗} = *C*^{∗}*V*^{∗}/*ρ* where *V*^{∗} = (*L*^{∗})^{3} is the volume unit and *ρ* the number density of tDPD particle. Here, the divisor *ρ* means that the unit volume is divided equally into *ρ* particles. Then, a concentration of 1.0 nM is scaled into the concentration unit of the tDPD system.

For example, in Section III C, species I (fibrinogen) has initial concentration of *C*_{I,0} = 7.0 × 10^{3} nM = 7.0 × 10^{−3} mol/m^{3}, the volume unit is *V*^{∗} = (*L*^{∗})^{3} = 1.0 × 10^{−18} m^{3}, and the number density of tDPD particle is *ρ* = 4. Thus, the moles of species I carried by each tDPD particle is *C*_{I,0} ⋅ *V*^{∗}/*ρ* = 1.75 × 10^{−21} mol, which should be scaled by *n*^{∗} = *C*^{∗} ⋅ *V*^{∗}/*ρ* = 1.0 nM ⋅ *V*^{∗}/*ρ*. Therefore, *C*_{I,0} = 7.0 × 10^{3} nM is mapped into tDPD unit as *C _{I}* = 7.0 × 10

^{3}.