key: cord-0784314-8scg3uer
authors: Yano, Ryosuke
title: Kinetic Modeling of Local Epidemic Spread and Its Simulation
date: 2017-03-16
journal: J Sci Comput
DOI: 10.1007/s10915-017-0408-9
sha: ed9f12e73f35ff02ca404add750abe36f592e8f9
doc_id: 784314
cord_uid: 8scg3uer

The local epidemic spread in physical space is modeled using the kinetic equation. In particular, the infection occurs via the binary interaction between the uninfected and infected individuals. Then, the local epidemic spread can be modeled on the basis of the stochastic Boltzmann type equation. In this paper, the normalized virus titer inside the infected human body is defined as the function of the elapsed time, which is measured from the infection time. Consequently, the probability of the infection at the binary human interaction increases, as the normalized virus titer inside the human body increases, whereas the normalized virus titer inside the infected human body decreases, after the normalized virus titer reaches to its maximum value, namely, unity, in the characteristic time. Numerical results indicate that the propagation speed of the boundary between the infected and uninfected domains depends on such a characteristic time, strongly, when the Knudsen number and temperature are fixed. Such a dependency of the propagation speed of the boundary between the infected and uninfected domains on the characteristic time can be described by the Fisher–Kolmogorov–Petrovsky–Piscounov equation which is introduced from the stochastic Boltzmann type equation. Finally, we consider three types of the human behavior as plausible actions to the local epidemic spread.

The epidemic spread sometimes does serious harms to humankind on the earth. In Tokio Marine and Nichido Risk Consulting Co., Ltd, 5-1, Otemachi, 1-Chome, Chiyoda-ku, Tokyo, Japan and do harms to people, previously. Thus, the quantitative prediction of the epidemic spread is markedly significant to suppress harms by the virus at minimum. The most classical mathematical-model to demonstrate the epidemic spread was proposed by Kermack and McKendrick [5] as SIR model (S: susceptible, I: infected, R: removed/recovered). In SIR model, numbers of individuals in statuses S, I, and R are temporally evolved. Similarly, the extended model of SIR model, such as SIRS with the birth and death [6] , SIRS model [7] and SIR model with other counterparts [8] , were proposed. These SIR type models are simultaneous ordinary derivative equations in terms of the time (t ∈ R + ). Thus, the epidemic spread in the physical space is not considered in SIR type models. To demonstrate the epidemic spread in physical space, the epidemic dynamics in the phase space (t, x) (in R + × X 3 , X 3 ⊆ R 3 ) must be considered. Cardy and Grassberger [9] considered the epidemic spread on network in the physical space using SIR model and found the similarity between the epidemic spread and percolation on the network. In recent studies by Barrat and his coworkers [10] [11] [12] [13] , the epidemic spread on the basis of the complex network theory has been investigated. Additionally, the lattice gas cellular automaton analysis of the epidemic spread was proposed by Fuks and Lawniczak [14] . For example, the epidemic spread in Japan is calculated using 69 nodes (major cities), where the number of inhabitants is large, adequately. Two nodes are connected, when the distance between two cities is shorter than the critical value. The individuals move from one node to other nodes, when they are connected each other. Figure 1 shows the fraction of the infected individuals (number of infected individuals at a node/total number of individuals at a node) at t = 0 and 10, in which t is the normalized time. At t = 0, the infected individuals exist in Okayama city. Readers can easily understand that the number of nodes with high fraction of infected individuals increases, as the time passes from the frame at t = 10. Theses studies of the epidemic spread in framework of the network theory or lattice gas cellar automaton are useful for understanding of the epidemic spread from a point (cell) to a point (cell) without tracing motions of infected individuals. Generally, the congested areas and depopulated areas sparsely exist in the country. Then, the epidemic spread from one node to the other node is reasonable approach. Meanwhile, we cannot explain the "local" epidemic spread inside a node using SIR model. Therefore, we must consider the local epidemic spread among individuals in the multi-agents system by tracing motions of individuals in order to investigate the local epidemic spread inside a node (city, village, county etc.).

Remembering that the harmful virus is transferred from one individual to other individuals, the simulation of the local epidemic spread can be performed by tracing motion of individuals, who are infected or uninfected. Provided that infections from infected individuals to uninfected individuals occur via binary interactions, exclusively, binary interactions occur, when the distance between two individuals is less than the two times of the effective diameter of an individuals, in which the effective domain of an individual is expressed with the volume of a sphere with an effective diameter (d). Then, d must be determined to reflect a droplet infection together with the contact infection. Such a long range interaction can be expressed with an intermolecular potential [15] . After all, the local epidemic spread on the basis of such motions of individuals can be formulated by the "stochastic" Boltzmann type equation. The kinetic equation such as the stochastic Boltzmann type equation is certainly useful for a calculation of a large number of individuals, which is hard to trace of motion of each of individuals [16] . Here, readers remind that the phrase "stochastic" means that a small number of individuals also can be discussed in the framework of the kinetic theory [17] .

A serious problem involved with SIR model is that transition rates such as S → I and I → R cannot be determined in the framework of the macroscopic scale. Provided that the infection status is expressed with the virus titer, we expect that a series of the infection process (uninfected → infected → recovery) can be understood with only one parameter, namely, the virus titer. A novelty of this paper is a kinetic modeling and computation of the epidemic spread using one parameter, namely, virus titer, whereas SIR type models distinguish individuals using three parameters (S, I and R). Readers readily understand SIR model allows the recovery of individuals, whose infection-interval is markedly shorter than the usual infection-term, because the transition of I → R can occur when a chosen random number (∈ [0, 1]) is smaller than the transition probability. Of course, an immunity is also another significant parameter to express the susceptibility of an individual toward the virus, which characterizes the probability of the infection of an individual. We, however, assume that the immunity is uniform for all the individuals to simplify our discussion. The virus tier, which is excited from zero by the binary interaction with the infected individual with some probability, is the state variable, which depends on the elapsed time after the infection. Then, the probability of the infection at the binary interaction is determined using the virus titer inside the infected individual. In particular, the normalized virus titer ( p ∈ [0, 1]) is defined as a function of the elapsed time (τ ), which is measured from the start time of the infection. Such a function must be determined to fit to real datum of the virus titer versus the elapsed time, as discussed by Beccam et al. [18] , whereas the normalized virus titer is defined as a Gaussian function of the elapsed time to simplify our discussion. Therefore, the realistic function of the normalized virus titer is set as our future study.

Finally, the distribution function in the stochastic Boltzmann type equation is defined using variables such as t, x, v (v ∈ V 3 : V 3 ⊆ R 3 : velocity space) and p (τ ). The primary aim of this paper is the comprehension of the characteristics of the local epidemic spread at the initial stage, which changes in accordance with variations of parameters in the stochastic Boltzmann type equation. In particular, we consider the local epidemic spread at the initial stage in the framework of the kinetic theory. Therefore, we never postulate the epidemic spread via vehicles such as cars, trains or airplanes, which were discussed by Colizza et al. [12] . Hence, the time scale of the evolution of the virus titer and time scale of motion of individuals are set unspecified parameters, because we never specify the species of the virus. Additionally, we never mention to multi-epidemic-sources and restrict ourselves to single epidemic-source (domain) to simplify our discussion. In numerical analyses of the stochastic Boltzmann type equation, the direct simulation Monte Carlo (DSMC) method [19] on the basis of the null time counter (NTC) method by Bird [19] is applied. The DSMC method to solve the stochastic Boltzmann equation, which demonstrates the epidemic spread on the basis of the virus titer, is surely new approach to analyses of the epidemic spread, whereas the numerical algorithm in the DSMC is similar to conventional one [19] which has been used to solve the Boltzmann equation for reactive gases, because the Monte Carlo judgment of the transition of the uninfected individual to the infected individual is similar to the Monte Carlo judgment of the chemical reaction via a binary collision in reactive gases [19] . Here, readers must remind that we solve the "stochastic" Boltzmann type equation, in which the phrase "stochastic" means that Klimontovich's distribution function [20] is considered in order to investigate effects of thermal fluctuations on the local epidemic spread. We introduce the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation [21] from the stochastic Boltzmann type equation in order to calculate the propagation speed of the boundary between the infected and uninfected domains. Finally, the local epidemic spread via three types of extended human behaviors (stochastic lethal/quarantine without interactions, hypokinetia, and refugee-quarantine with interactions) in accordance with the local epidemic spread is investigated. One advantage of the kinetic model based on the virus titer over SIR model is that we can construct a flexible mathematical-model to reflect scenarios of actions of individuals in accordance with the virus titer of the infected individual (i.e, quarantine, refugee, hypokinesia etc.). In this paper, we primarily focus on the local epidemic spread owing to the thermal motion of individuals which does not postulate the intended velocity of pedestrians, because our aim of the study is an investigation of the local epidemic spread on the basis of the virus titer coupled with thermal motions of individuals in the framework of the kinetic theory. Thus, further advanced interaction models beyond the isotropic scattering in the binary collision between two pedestrians, which were proposed by Helbing [15] or Gipps and Marksjö [22] , are not discussed in this paper, because the author considers that the investigation of the local epidemic spread, which includes the detailed dynamics of pedestrians, seems to be out of scope of this study. We, however, mention to the local epidemic spread, in which intended velocities of pedestrians are fixed to their constant values in "Appendix" in order to prove that the proposed kinetic model can be extended to more realistic model of the local epidemic spread in our future study. Finally, the Bhatnagar-Gross-Krook (BGK) like term is added to the stochastic Boltzmann equation to describe the relaxation of velocities of pedestrians toward their intended velocities. This paper is organized as follows. In Sect. 2.1, the stochastic Boltzmann type equation for the local epidemic spread via the binary collision between the infected and uninfected individuals is formulated. In Sect. 2.2, the propagation speed of the boundary between the infected and uninfected domains is calculated from the FKPP equation, which is introduced from the stochastic Boltzmann type equation. In Sect. 2.3, DSMC analyses of the local epidemic spread are performed to investigate effects of parameters, which are used in the stochastic Boltzmann type equation, on the local epidemic spread. In Sect. 3, three types of the human behavior in accordance with the local epidemic spread are modeled and characteristics of three types of the human behavior in accordance with the local epidemic spread are investigated, numerically. Finally, we make concluding remarks in Sect. 4.

Numerical Analysis 

The normalized virus titer inside the human body is expressed with p, in which p = 0 corresponds to the uninfected state and 0 < p corresponds to the infected state. Here, we remind that p := p (τ ) ∈ R + is a function of the elapsed time (τ ), once an individual is infected by the interaction with the infected individual at t = t p (t, t p ∈ R + ) such as

where τ := t − t p ∈ R + is the elapsed time, after the virus intrudes into the human body. From Eq. (1), p(0) > 0 and p (τ ) ∈ (0, 1] is obtained. τ m corresponds to the characteristic time, when the normalized virus titer reaches to the maximum value and T c is the variance of p (τ ).

The time derivative of p (τ ) is obtained as

The motion of the individual is characterized by the translation in physical space X 3 ⊆ R 3 and binary elastic collision (interaction) with another individual, which yields the change of the direction of his or her motion. Therefore, such a random change of the direction of motion via the binary elastic collision between two individuals seems to be implausible in our real society, because each of individuals tends to move toward his or her destination in accordance with his/her intended velocity. Meanwhile, we consider that such a random change of the direction of motion seems to be plausible, when children move in the classroom or campus of the school, as they like. The virus inside the body of the infected individual intrudes into the body of the colliding partner in accordance with the following probability:

where W ∈ [0, 1] is the white noise. Equation (3) indicates that the probability of the virus-infection increases, as p (τ ) increases. Thus, the probability of the virus-infection is low at the incubation period (τ τ m ) or complete recovery period (τ m τ ). The Klimontovich's distribution function [20] is defined by

where the subscript "i" is an index of the i-th individual, N c is the total number of individuals in X 3 , V c := X 3 v is the infinitesimal volume in x ∈ X 3 v (X 3 v ⊆ X 3 ) and V 3 ⊆ R 3 : velocity space.

From above discussion, the stochastic Boltzmann type equation, which demonstrates the infection of the virus owing to the binary elastic collision between the infected and uninfected individuals, is formulated as:

I 11 10 + I 10 10 + I 11 01 + I 01 01 + I 00 00 + I 11 11 Ẽ ξ σ sin χdχd d 3 v 1 dp 1 ,

where A ∈ R + has the unit [m/s] and fixed to 1 [m/s] in this paper, δ i j is Kronecker's delta function, and an index "s" in f (s) indicates the status of the infection, namely, s = 0 (non-infection) and s = 1 (infection). In Eq. (5), we assume that all the individuals have same mass (m) in order to simplify our discussion. An index "s" is, however, not a dependent parameter, because the status of the infection can be judged by the normalized virus titer, namely, p(τ ) = 0 (non-infection) and p(τ ) > 0 (infection) from Eq. (1). The subscription "1" indicates the collisional partner. In the right hand side of Eq. (5),Ẽ is defined as

whereṽ andṽ 1 are normalized velocities by C ∞ := √ 2Rθ ∞ such asṽ := v/C ∞ and v := v 1 /C ∞ (R: gas constant, θ ∞ : temperature).

From Eq. (5), the collision frequency depends on the total collisional energy, because an individual with the higher kinetic energy tends to access to another individual the more frequently, whereas the collision frequency in the Boltzmann equation depends on the relative velocity (|v − v 1 |). Here, we are free from the geometric collisional mechanics, which is postulated in the Boltzmann equation. Additionally, ξ ∈ [0, 1] is a parameter, which expresses the dependency of the collision frequency on the total collisional energy. σ , χ ∈ [0, π] and ∈ [0, 2π] in the right hand side of Eq. (5) are the differential cross section (σ = d 2 /4), deflection angle and scattering angle, respectively. Indeed, above isotropic scattering in the binary collision must be modified in order to demonstrate the characteristics of pedestrians, as discussed by Helbing [15] .

From the left hand side of Eq. (5), f ∈ C 1 X 3 ∧ f ∈ C 1 V 3 is postulated, although Klimontovich's distribution function does not always satisfy the smoothness and continuity of f in both X 3 and V 3 , whereas ψ ∈ C 1 V 3 is always satisfied by d p ψ( p) = −c p (∵ c p < +∞) from Eq. (2).

Multiplying 1, p,..., p n (n ∈ N) by both sides of Eq. (5) and integrating over V 3 × P, we obtain

where G p n := V 3 ×P p n Q ( f, f 1 ) d 3 vdp is obtained as

and ρ p ,p n , δu p n are defined as

From Eqs. (7) and (8), the propagation velocity of the virus titer moment (ρ ppn ) is δu p n . The number density of the infected individual ρ p increases via the collisional term G p 0 . Therefore, the propagation of the infected domain depends on G p 0 together with the term ∇ x · ρ pp0 δu p 0 .

In the left hand side of Eq. (7), terms c p ρ ppn (1 ≤ n) satisfy

Equation (10) indicates that the virus titer moment increases at 0 ≤ τ < τ m and deceases at τ m ≤ τ owing to the time evolution of the virus in Eq. (1)

.

, the propagation velocity of the virus titer moment never depends on the distribution of p. We, however, conjecture that an individual with the larger kinetic energy can be a seed of the infection at the uninfected domain the more frequently. In later discussion on numerical results in Sect. 2.3, we confirm that δu p 0 2 is surely larger than |δu| 2 at the boundary between the infected and the uninfected domains.

Besides with virus titer moments in Eq. (9), we define some moments such as

In Sect. 2.1, we formulated general form of moment equations for the infected individuals from Eq. (5). On the other hand, moment equations for the uninfected individuals are formulated for individuals with p = 0, namely, p ∈ P − P.

Multiplying 1 by both sides of Eq.

(1) and integrating over P − P × V 3 , we obtain

where G un := (P −P)×V 3 I 11 01 + I 01 01 + I 00 00 dvd p = −G p 0 is used. From Eqs. (7) and (12), we readily confirm the mass conservation, namely, ∂ t ρ +∇ x · ρ pp0 δu p 0 + ρ un δu un = 0. The essence of the diffusion flux (J) is an expression of the separated velocity betweenp 0 δu p 0 and δu p un with δu.

Using J, Eqs. (7) and (12) can be rewritten as

where J is defined using δu := ρ −1 ρ p 0 δu p 0 + ρ un δu un in Eq. (11) as

In the framework of the Navier-Stokes-Fourier approximation, the diffusion flux is approximated as [23]

where η p 0 := ρ p 0 ρ −1 is the mass fraction of the number of infected individuals to the total number of individuals at x, D is the diffusion coefficient and D p is the pressure diffusion coefficient, m is a mass of the individual and p := ρ Rθ is the static pressure. G p 0 can be written as

where K 01 is the transition rate by which the uninfected individuals change to infected individuals.

Setting δu = 0, m = 1, D p = 0 and D =const. in X 3 , we obtain the following mass transport equation from Eqs. (13) , (15) and (16):

Equation (17) coincides with the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation [24] , when ρ = 1, as follows

Consequently, the epidemic spread has an analogy to the dynamics described by the FKPP equation, when forgoing assumptions used to obtain Eq. (18) from Eq. (5) is validated. On the other hand, Succi and his coworkers analyzed the FKPP equation with nonzero flow-flux, namely, δu = 0 using the lattice Boltzmann method (LBM) [25, 26] . The application of the LBM to the FKPP equation, which demonstrates the epidemic spread, is interesting issue in our future study, definitely. The reaction-diffusion equation can be surely applied to the conventional SIR model using similar procedures mentioned above such as

where η S ∈ R + (fraction of susceptible individuals), η I ∈ R + (fraction of infected individuals), η R ∈ R + (fraction of recovered individuals), and K i→ j ∈ R + is the transition rate from status "i" to " j". Figure 2 shows the time evolution of η I , which is obtained by solving Eq. (19) using D = 1,K S→I = 0.25 andK I→R = 0.1, where all the physical quantities are normalized. Thus, the normalized quantities are denoted by˜. x ∈ X ⊆ X 2 (x ∈ X||x| ≤ 100 ∧ |ỹ| ≤ 100) is considered as the numerical domain with the periodic boundary. The finite differential scheme is used to solve Eq. (19) , whereas Succi [26] indicated that FKPP equation with the advection can be solved with better accuracy using the LBM rather than the finite differential scheme. Actually, we consider the directed flow of individuals in the process of the epidemic spread, as discussed in Sect. 3 and appendix. As an initial condition, η S = 0.9, η I = 0.1 and η R = 0 in the domain X I (x ∈ X I | x 2 +ỹ 2 ≤ 50), whereas η S = 1 and η I = η R = 0 in the domain X − X I . The infected domain expands, concentrically, ast increases. Most of infected individuals at X I andt = 0 change to recovered individuals at X I andt = 30. Finally, most of infected individuals change to recovered individuals att = 200.

The propagation speed of a traveling wave [24] which is obtained by Eq. (18) can be calculated by following procedures, when f in X 3 is spherically symmetric, [27] 1. We setφ := η p 0 andψ := dη p 0 /dz using the definition z := r − ct (r := ||x|| L 2 ). 2. Equation (18) can be rewritten usingφ andψ as 

whereφ := dφ/dz andψ := dψ/dz. 4. Two eigenvalues in Eq. (21) are obtained as

5. We assume that the boundary between two domains X 3 p ( x ∈ X 3 p |η p 0 (x) ∈ (0, 1] ) and X 3 − X 3 p , namely, infected and uninfected domains, is expressed with φ ,ψ = (0, 0). The propagation speed of the boundary between such two domains must satisfy the following relation to yield the real value of λ ± at φ ,ψ = (0, 0):

From item 5, the minimum value of the propagation speed of the boundary φ ,ψ = (0, 0)

is equal to 2 √ DK 01 . In above numerical result of the reaction-diffusion equation of SIR model in Eq. (19) , the propagation speed of the boundary between the infected and uninfected domains is defined by 2 DK S→I = 1. Figure 3 shows the time evolution of η I . The propagation speed of the boundary certainly coincides with unity, as shown in Fig. 3 . Provided that D is independent of τ m and T c in Eq. (1), the minimum propagation speed of the boundary between the infected and uninfected domains, namely, 2 √ DK 01 depends on τ m and T c via K 01 owing to Eq. (8) . In this paper, we never refer to concrete form of K 01 owing to mathematical complexities, whereas later numerical results show that the propagation speed of the boundary domain strongly depends on τ m , when T c is fixed. From Eq. (8), we make a rough estimation of K 01 as K 01 ∼ K p (τ ) (K ∈ R + ). Figure 4 shows

p|τ m =3/4 / p|τ m =1 decreases from 8.91 atτ = 0 in the case ofτ m = 3/4, asτ increases. Then, we can understand that the propagation speed of the boundary domain between the infected and uninfected domains is markedly sensitive toτ m .

All the physical quantities are normalized as 

Here, the distribution function in x 3 -axis is assumed to be uniform. Of course, (v, x) ∈ V 2 ×X 2 might be more realistic to demonstrate local epidemic spread between the infected and uninfected individuals than (v, x) ∈ V 3 ×X 3 , when we exclude the interaction between two individuals inside three dimensional facilities (i.e., binary collision at the staircase.), whereas we consider (v, x) ∈ V 3 × X 3 , because the assumption of the uniformity along x 3 -axis does not change the qualitative tendency of the local epidemic spread in (x 1 , x 2 ) space, markedly, as the two dimensional flow has been considered using not two dimensional collision between two discs but three dimensional collision between two spheres [28] . The reduction of the dimension will be easily done in accordance with the situation to be considered. Here, we setx := X ,Ỹ ,Z and calculation domain as X ,Ỹ ,Z ∈ X ⊆ X 3 , whereas we set δũ := (δũ, δṽ, δw) and δũ p n := δũ p n , δṽ p n , δw p n . As an initial condition, we set T c /t ∞ = 0.05,ρ = 1,ũ ∞ = 0, and Kn = ( √ 2πd 2 ρ ∞ /m L ∞ ) −1 = 1.0 (m: common mass of individuals). N c = 100 sample particles per a lattice are set, in which 61 × 61 lattices equally spaced inside the calculation domain, X, {x,ỹ ∈ X|x ∈ [0, 1] ,ỹ ∈ [0, 1]}. As the boundary condition onX = 0, 1 and Y = 0, 1, the periodic boundary condition is used in Tests A-1, 2, 3, 4, 5, B and C, whereas the specular boundary condition is used in Test A-6. Att = 0, all the individuals inside the

In other words, all the individuals have p (τ = 0) > 0 owing to Eq. (1) inside the domain U att = 0, whereas individuals inside the domain X − U, who have p = 0 att = 0, are uninfected. The schematic of the calculation domain (X) is shown in Fig. 5 (Table 1 ). Figure 6 shows the time evolution ofp in Test A-1. The infected domain spreads at 0 ≤t ≤ 1.0, concentrically, where the center of the concentric spread is, certainly, X ,Ỹ = (0.5, 0.5). Atτ m = 1.0 in Eq. (1), p of the individual, who is infected att = 0, reaches to unity atτ = 1. Consequently,p increases toward X ,Ỹ = (0.5, 0.5) and reaches to 0.895 at X ,Ỹ = (0.5, 0.5) att = 1, where we remind thatp is not equal to unity at t = 1, because individuals, who are infected att > 0 and have p (τ ) < 1 owing toτ < 1, (1) must be damped to zero, rapidly, before the interaction with the uninfected individuals. Such a situation can be obtained by using p (τ c ) 1 (τ c : inverse of the averaged collision frequency), whenτ m <τ c . Next, we investigate the effect ofτ m on the local epidemic spread by considering Test B. Figure 7 shows the time evolution ofp in Test B. The infected domain spreads, concentrically in a similar way to Test A. The propagation speed ofp in Test B is higher than that in Test A-1. From Eq. (7), the propagation of the infected domain depends on the spatial propagation via ∇ x ρ 0p0 δu p 0 and generation of ρ 0p0 via G p 0 . We remind thatp is related (1) rather than the increase in δu p 0 , because our numerical results indicate that δũ p 0 t at the boundary between the infected and uninfected domains in Test A-1 is similar to that in Test B. Such a similarity between δũ p 0 t at the boundary between the infected and uninfected domains in Test A-1 and that in Test B indicates that D in Test A-1 is similar to that in Test B. In short, the difference between the propagation speed of the boundary domain between the infected and uninfected domains in Test A-1 and that in Test B is caused by the relation K 01 | Test A-1 K 01 | Test B as a result of Eq. (23). Next, we consider on effects of velocity fluctuations via thermal fluctuations by comparing ofp, which are obtained in Tests A-1, 2 and 3. In the DSMC method to solve the stochastic Boltzmann type equation, the velocity fluctuations (δu) depend on the number of sample particles in a lattice (N c ), which is averaged during some time interval, such that [29] 

whereθ is the temperature, which is averaged during some time interval. Figure 8 shows snapshots ofp in Tests A-1, A-2 and A-3 att = 1.0 (upper frames) and t = 1.536 (lower frames). Figure 8 shows that there are not marked differences between contours ofp in Test A-1 and those in Test A-2 att = 1.0 and 1.536, whereas there are marked differences between contours ofp in Test A-1 and those in Test A-3 att = 1.0 and 1.536, because the infected domain never spreads concentrically in Test A-3. From Eq. (24), we can conjecture that δũ 2 in Test A-1 is a half of that in Test A-2 and 1/10 times of that in Test A-3, becauseθ is fluctuating around 5 at 0 ≤t in Tests A-1, 2 and 3. The dependency of the propagation speed of the boundary between the infected and uninfected domains on δũ is ignorable, because there are no marked differences among the area of the infected domain in Test A-1 and those in Tests A-2 and A-3 att = 1 and 1.536, as shown in Fig. 8 . Consequently, the propagation speed of the boundary between the infected and uninfected domains strongly depends on the collisional term G p 0 in Eq. (7), in other words, K 01 in Eq. (23) . Next, we investigatef (τ ) := V 3f t,x, p (τ ) ,ṽ dv 3 . We set τ = 0.02 in the range of 0 ≤τ ≤ 2.f (τ ) is the set of simple functionsf ( τ × i) (i ∈ Z + ) on Lebesgue measurẽ τ ∈ i × τ , (i + 1) × τ ). Then, 0 ≤τ ≤ τ includes both uninfected individuals and infected individuals, whose elapsed time is in the range of 0 ≤τ ≤ τ . Figure 10 shows in Test C occur more frequently than that does in Test B, because the propagation speed of the boundary between the infected and uninfected domains in Test C is lower than that in Test B. Consequently, the range ofτ , which yieldsf (τ ) = 0, namely, sum of Lebesgue measure, which yieldsf (τ ) = 0, becomes the longer, asτ m becomes the longer. In summary, the ratio of the number of individuals, who are infected at other lattices (domains), to the total number of infected individuals at the specific lattice (domain) decreases, as τ m decreases, because the local epidemic spread occurs inside the lattice rapidly via one or a few of invasions of the infected individuals from neighboring lattices, before some infected individuals reach to neighboring lattice, because the virus titer increases rapidly inside the lattice owing to small τ m from Eq. (2). Figure 12 shows the sample particles (individuals), which are calculated using the DSMC method, in Test C (att = 0.6) and Test A-1 (att = 0.88). A sphere corresponds to one individual, where the diameter of a sphere is expressed with enlarged Fig. 12 . Consequently, we confirm that a smallτ m yields the rapid local epidemic spread via the binary collision between the infected and uninfected individuals, once a infected individual invades to the outside of the boundary between the infected and uninfected domains. Thus, the infected individual with the large velocity becomes a seed of the local epidemic spread, which will be confirmed in the discussion of δu p 0 at the boundary between the infected and uninfected domains in Fig. 13 . Consequently, such a higher velocity of the infected individuals at the boundary between the infected and uninfected domains attributes to the finite diffusion flux from Eq. (14) at the boundary between the infected and uninfected domains.

Toward further investigation of the propagation speed of the boundary between the infected and uninfected domains, we calculate the propagation speed of the traveling wave front in Tests A-1, B and C. The rough estimation of the propagation speed of the traveling wave front can be performed from Figs. 10 and 11 by calculating the arrival time of the wave front at point (B) X ,Ỹ = (0.5, 0.25) such asṽ * |τ m =1 = 1.25,ṽ * |τ m =0.75 = 2.5 andṽ * |τ m =0.5 = 12.5. As a result, we obtainṽ * |τ m =0.75 /ṽ * |τ m =1 = 2 andṽ * |τ m =0.5 /ṽ * |τ m =1 = 20. Orders of such two ratios of the traveling wave front are similar to those calculated using Eq. (23), as shown in Fig. 4 . Finally, the propagation speed of the traveling wave front depends onτ m strongly, as predicted by discussions in Sect. 2.2.

Next, we investigate effects of the dependency of the local epidemic spread on ξ in the collision term in Eq. (5) andθ , respectively. Figure 13 shows snapshots of the vector field (ũ p 0 ,ṽ p 0 ) in Tests A-1 (left frame) and A-4 (middle frame) att = 1.0 and A-5 (right frame) att = 0.72, whereũ p 0 := (ũ p 0 ,ṽ p 0 ,w p 0 ). The propagation speed of the boundary between the infected and uninfected domains in Test A-4 is slightly lower than that in Test A-1, whereas the propagation speed of the boundary between the infected and uninfected domains in Test A-5 is markedly higher than that in Test A-1. As a result, the increase in ξ yields the slight increase in the propagation speed of the boundary between the infected and uninfected domains owing to the increase in the collision frequency, because the majorant collision frequency [30] is determined using 1 <Ẽ max in the DSMC calculation. Additionally, the increase inθ yields the increase in δu p 0 at the boundary between the infected and uninfected domains, which leads to more rapid propagation of the boundary between the infected and uninfected domains in Test A-5 than that in Test A-1. ig. 14 Snapshots ofp att = 1.0 in Tests A-6 (upper-left frame) and A-1 (upper-right frame), and at t = 1.536 in Tests A-6 (lower-left frame) and A-1 (lower-right frame) uals with the large velocity, who determine the propagation speed of the boundary between the infected and uninfected domains, increases, asθ increases. Actually, we can confirm that the propagation speed of the boundary between the infected and uninfected domains depends on not δũ but δũ p 0 at the boundary between the infected and uninfected domains from Figs. 7 and 13. The increase in δu p 0 at the boundary between the infected and uninfected domains owing to the increase inθ can be explained by the fact that D increases in accordance with the increase inθ due toD ∝θ ξ/2 in Eq. (14) . Meanwhile, K 01 also increases in accordance with the increase inθ , when ξ > 0. Thus, it is still open problem whether the increase in the propagation speed of the boundary between the infected and uninfected domains via the increase inθ is caused by the increase in D or K 01 in Eq. (23) . Meanwhile, the description of the relation between δũ and δũ p 0 at the boundary between the infected and uninfected domains is set as our future study. In Fig. 13 , we find that δu p 0 at the boundary between the infected and uninfected domains tend to be lager than δu p 0 inside the infected domain. Such a tendency indicates that the virus tier is correlated to the velocity of the individual such as f ( p, v) C ( p, v) ) is the correlation function).

Finally, we investigate effects of the boundary condition on the local epidemic spread. Figure 14 shows snapshots ofp att = 1.0 in Tests A-6 (upper-left frame) and A-1 (upperright frame), and att = 1.536 in Tests A-6 (lower-left frame) and A-1 (lower-right frame). There are no marked differences betweenp in Test A-6 andp in Test A-1 att = 1.0 and 1.536. Consequently, effects of the boundary condition on the local epidemic spread are ignorable.

From forgoing discussions on numerical results, we obtain some significant conclusions on the local epidemic spread as follows:

(i) The propagation speed of the boundary between the infected and uninfected domains is dominated by bothτ m andθ , exclusively. (ii) The boundary between the infected and uninfected domains spreads rapidly, whenτ m is small. Then, the infection with smallτ m can spread by a few invasions of infected individuals to the uninfected domain, when the interaction between the infected and uninfected individuals is allowable. Consequently, we must pay attention to the local epidemic spread with the smallτ m by suppressing movements of infected individuals not to invade to the uninfected domain. Of course, the local epidemic spread can be accelerated by the transport of the infected individuals to the uninfected domain using the vehicle, whenτ m is small. 1 (iii) The infected individual with the large velocity can invade to the uninfected domain, more frequently than those with the small velocity. The higher velocity of the infected individuals than that of the uninfected individuals at the boundary between the infected and uninfected domains yields the finite diffusion flux, which is expressed by Eq. (14) . Then, the infected individual with the large velocity can be a seed of the infection in the uninfected domain, more frequently than those with the small velocity. Consequently, the virus titer p is correlated to |v| in f (t, p, v, x) . The diffusion flux resulted from the difference between such two velocities at the boundary between the infected and uninfected domains enables us to demonstrate the dynamics of the epidemic spread using the FKPP equation in Eq. (18) . Actually, the dependency of the propagation speed of the traveling wave front (i.e., boundary between the infected and uninfected domains) onτ m is described by Eq. (23) with some good accuracies.

In above discussion, the probability of the infection depends on p > 0 (or p 1 > 0) andẼ, as shown in Eq. (5) . Here, we consider the extension of the infection process of the virus by proposing three models. The first model (model A) is that the infected individuals with the normalized virus titer p are excluded from the calculation domain with the rate μp. In other words, infected individuals pass away or quarantined with the rate μp.

In the model A, the infected individuals with p = 1 lead to death or being quarantined with the maximum probability.

Consequently, Eq. (5) is modified as

where μ ∈ R + corresponds to the force of mortality or quarantine. Provided that μ is equal to the inverse of the time interval 1/ t, the infected individuals with p = 1 always lead to death or being quarantined. In this paper, we restrict ourselves to the case of μ = 1/ t to investigate effects of μ on the local epidemic spread.

As the time interval, four types of the time interval are considered such as t = 2 × 10 −3 (Test D-1), 4 × 10 −3 (Test D-2), 5 × 10 −4 (Test D-3) and 1 × 10 −4 (Test D-4). Other initial conditions in Tests D-1, 2, 3 and 4 are quite sames as those in Test A-1. Figure 15 shows time evolutions ofρ (top six frames),θ (middle six frames) andp (bottom six frames) in Test D-1. The infected individuals in U pass away (or are quarantined), as shown in contours ofρ att = 0.0002, 0.8 and 1.6. Meanwhile, the number of individuals, who pass away or are quarantined, is large alongX = 0 andỸ = 0, as shown in contours ofρ att = 2.4 and 3.2. Finally, the infected individuals have disappeared att = 8.8. Therefore, Model A allows the survival of individuals, who are not infected, by arranging the force of mortality or quarantine, namely, μ, as discussed later.θ fluctuates around 5 att = 0.002, whereas there emerges hot spots around X ,Ỹ = (0.5, 0.5) att = 0.8, 1.6 and 2.4. From contours ofρ andθ , we can confirm thatθ is high in the domain, whereρ is small, att = 0.8, 1.6, 2.4 and 3.2. The averaged temperature in the calculation domain, namely, (1/V ) Xθ d 3 x (V : volume of the calculation domain), temporally decreases fromt = 2.4 to 3.2. Such a decrease inθ indicates that individuals with the higher kinetic energy tend to pass away or being quarantined the more frequently, because the infection rate via the binary collision increases in accordance withẼ ξ , as shown in Eq. (5) . Figure 15 shows that infected domain spreads concentrically att = 0.8, as shown in the contour ofp, whereasp is large locally att = 1.6, 2.4 and 3.2. We confirm thatp is large in the domain, whereρ is small, owing to the higher probability of the death or quarantine in accordance with the larger p in Model A. The infected individuals have disappeared att = 8.8, as shown byp = 0 in X at t = 8.8.

Finally, we investigate effects of μ on the total number of uninfected individuals, namely, Figure 16 showsρ t versust in Tests D-1, 2, 3 and 4, which are obtained using Model A. Figure 16 shows the result: These results indicate that limt →+∞ρt t increases, as μ increases. Thus, the more rapid death or quarantine, which exhibits infected individuals to interact with others, yields the more uninfected individuals att → +∞. Meanwhile, such a rapid death or quarantine via the increase in μ shows the rapid decrease ofρ t in the early term (0 ≤t < 0.6), as shown in the enlarged view in Fig. 16 . Indeed, we can understand that the rapid quarantine of the infected individuals via the early detection of the virus is significant to avoid the local epidemic spread. In the second model (model B), the velocity of the infected individuals decreases in accordance with p such as

where ζ p is the damping rate. 

The total collisional energy (Ẽ) decreases owing to the decrease in the velocity of the infected individual in model B. Consequently, the infection rate via the binary collision decreases, as the damping rate of v, namely, pζ p , increases in accordance with the increase in p (τ ).

To investigate the characteristics of Model B, the local epidemic spread is calculated using ζ p = 10, whereas initial conditions are quite same as those in Test A-1. Figure 17 shows time evolutions ofρ (top six frames),θ (middle six frames) andp (bottom six frames), which are obtained using Model B. The infected domain spreads concentrically and the propagation speed of the boundary between the infected and uninfected domains at 0 ≤t ≤ 1, which is obtained using Model B, is similar to that in Test A-1, as shown in contours ofp in Figs. 5 and 17. The decrease in the kinetic energy of the infected individual yields the low temperature domain around X ,Ỹ = (0.5, 0.5) att = 0.768. Such a decrease inθ yields the decrease in the pressure, which attributes to the inflow of individuals, who exist in the high pressure domain, to the low pressure domain. Consequently,ρ temporally increases around X ,Ỹ = (0.5, 0.5) at 0.768 ≤t. Finally, the concentric contour ofp att = 1 changes to the rhombus contour ofp att = 1.92 owing to the periodic boundary condition. Indeed, Model B fails to delay the propagation of the boundary between the infected and uninfected domains, because the decrease in the kinetic energy of the infected individual yields the low pressure regime, which attracts uninfected individuals from the high pressure regime. Finally, in the third model, uninfected individuals try to go away from the gravitational center (x G ) of infected individuals in accordance with the repulsive force (F r ) and infected individuals are gathered to their gravitational center in accordance with the attractive force (F a ). Such a situation still allows an interaction between infected and uninfected individuals. Therefore, gathering of the infected individuals is not complete quarantine, which removes existences of infected individuals with some probability from X like Model A. The repulsive and attractive forces are defined as

From Eq. (28), we obtain

As a result of Eq. (29), Eq. (5) is modified as

Model C (quarantine with interactions and refugee model)

To investigate the characteristics of Model C, numerically, κ r = 10 and κ a = 100 are set. Other initial conditions are same as those in Test A-1. Figure 18 shows time evolutions ofρ (top six frames),θ (middle six frames) andp (bottom six frames), which are obtained using Model C with κ r = 10 and κ a = 100. Contours ofp in Fig. 18 indicates that the propagation speed of the boundary between the infected and uninfected domains, which is obtained using the Model C with κ r = 10 and κ r = 100, is lower than that in Test A-1 owing to refugees of uninfected individuals and swarming of infected individuals around X ,Ỹ = (0.5, 0.5).

The infected individuals try to go away from infected individuals att = 1.0, as confirmed by the vector fields of (ũ,ṽ) in Fig. 19 . Additionally, velocities of individuals around the boundary between the infected and uninfected domains is large att = 2.16 and 3.98, whereas individuals cluster toward X ,Ỹ = (0.5, 0.5) inside the boundary between the infected and uninfected domains. As a result,ρ around X ,Ỹ = (0.5, 0.5) temporally increases, whereasθ inside such a clustering domain temporally increases at 0.24 ≤t, as shown in Fig. 18 . Finally, all the individuals in X are infected in Model C, although the propagation speed of the boundary between the infected and uninfected domains in Model C is lower than that in Test A-1.

The stochastic Boltzmann type equation was formulated to demonstrate the local epidemic spread in the local domain at the initial stage by focusing on the time evolution of the virus titer inside the infected individual. The propagation speed of the boundary between the infected and uninfected domains strongly depends on the characteristic time of the virus titer, namely, τ m , whereas the individual with the higher velocity becomes a seed of the local epidemic Vector fields (ũ,ṽ) att = 1.0, 2.16 and 3.98, which are obtained using Model C with κ r = 10 and κ a = 100 spread at the boundary between the infected and uninfected domains, the more frequently, as τ m decreases. The analytical result of the propagation speed of the boundary between the infected and uninfected domains is obtained using the FKPP equation. The dependency of the propagation speed on τ m , which is analytically estimated using the FKPP equation, is numerically confirmed with some good accuracies. Additionally, the absolute value of the fluctuating flow velocity of infected individuals, namely, δ|u p 0 | , at the boundary between the infected and uninfected domains increase in accordance with the increase in the initial temperature. The increase in δ|u p 0 | yields the increase in the propagation speed of the boundary between the infected and uninfected domains. Such an increase in the propagation speed of the boundary between the infected and uninfected domains owing to the increase in the initial temperature can be explained by increases in D and K 01 (0 < ξ) in the framework of the FKPP equation, as shown in Eq. (23) . The ratio of the number of individuals, who are infected at other lattices (domains), to the total number of infected individuals at the specific lattice (domain) decreases, as τ m decreases, because the local epidemic spread occurs inside the lattice rapidly via one or a few of invasions of the infected individuals from neighboring lattices, before some infected individuals reach to neighboring lattice, because the virus titer increases rapidly inside the cell owing to small τ m . Meanwhile, the change of the velocity fluctuations via the change of the number of sample particles does not affect the propagation speed of the boundary between the infected and uninfected domains. The dependency of the collision frequency on the total kinetic energy of two interacting individuals indicates that the stronger dependency of the collision frequency on the total kinetic energy yields the faster propagation of the boundary between the infected and uninfected domains owing to the increase in K 01 (ξ > 0). Finally, three types of the human behavior in accordance with the local epidemic spread was considered. In Model A, the more rapid death or quarantine of infected individuals, which is done not to interact with other uninfected individuals, yields the more uninfected individuals under t → +∞. In model B, the hypokinesia of infected individuals yields the concentration of individuals around the center of the infected domain. Consequently, Model B fails to delay the local epidemic spread owing to the increase in the collision frequency between the infected and uninfected individuals around the center of the infected domain. In Model C, the quarantine of the infected individuals inside the specific area was expressed with the attractive force, whereas the refugee of the uninfected individuals from the infected individuals was expressed with the repulsive force. Indeed, Model C succeeded the delay of the local epidemic spread, whereas the allowance of interactions between the infected and uninfected individuals leads to the result that all the individuals in the calculation domain are infected. Various types of the combination of Models A, B and C together with the realistic dynamics of pedestrians in "Appendix" are interesting in our future study together with the comparison of numerical results, which are obtained using our kinetic model, with real datum of the local epidemic spread.

In this paper, we have considered the dynamics of pedestrians (individuals), who never have their intended velocities [15] . Additionally, we have assumed the isotropic scattering after the binary collision, which might be rough model to describe the characteristics of the dynamics of the pedestrians [15] . For example, the number of pedestrian, who go back to the direction, which has been walked, after the binary collision is presumably less than that of pedestrian, who go straight toward the intended direction after the binary collision. Therefore, we must set the weight on the small deflection angle in order to emphasize the grazing collision, which does not change intended velocities of colliding two pedestrian. Such stress on the small deflection angle is attained by changing a weighting parameter in the deflection angle in the Variable Soft Sphere (VSS) model [31] . Additionally, the kinetic equation beyond the stochastic Boltzmann equation (i.e., Enskog equation) must be considered to demonstrate packing effects via the high volume fraction in X 3 occupied by pedestrian [32] . Finally, the sophistication of the collisional process in the local epidemic spread, which considers the characteristics of pedestrian, will be set as our future study.

Here, we mention to the effect of the intended velocities, briefly. The velocity of the individual with the index-i relaxes to intended velocities such as

where v 0 i is the intended velocities and ν is the relaxation rate. From Eq. (31), we can formulate the stochastic Boltzmann equation with the relaxation to intended velocities such as

where ψ 0 (t, v, x, p, s) is defined as

The relaxation term in Eq. (32) is a form of the BGK model, whereas Eq. (32) indicates that the energy-momentum conservation violates. The author recommends readers to access to another type of the Boltzmann type equation proposed by Helbing [33] , which demonstrates the collisional dynamics between two pedestrians by setting the intended velocity as a dependent variable in the distribution function. In our numerical tests in Sects. 2.2 and 3, the initial distribution function is obtained using the Maxwell-Boltzmann distribution. Similarly, we consider the numerical condition att = 0, which is quite same as that in Test A-1. Additionally, we assume that the intended velocities of pedestrian are equal to those att = 0. Figure 20 showsp att = 1.34, which are obtained using Eq. (32), whenν = 0.25 (left frame), 0.5 (middle frame) and 1 (right frame), in whichν := νt ∞ . Figure 20 shows that the propagation speed of the boundary between the infected and uninfected domains decreases, asν increases. As a result, we find that the intended velocity is a key parameter to determine the propagation of the boundary between the infected and uninfected domains.

Next, we refer to the limiting case, in whichν is markedly large and the intended velocities of pedestrian are limited to two values, namely, ±ṽ 0 , 0, 0 . The number of pedestrians with the intended velocity ṽ 0 , 0, 0 is same as that of pedestrians with −ṽ 0 , 0, 0 . Of course, the intended velocity of each of pedestrians is fixed to its constant value, namely, ±ṽ 0 , during the calculation. Figure 21 shows time evolutions ofp, whenṽ 0 = 1 andν = 100 in Eq. (32) , and initial numerical condition is quite same as that in Test A-1.

Three characteristics domains, S 1 , S 2 and S 3 , which are defined as are obtained att = 1.p in S 1 and S 2 is larger than that in S 3 , whereasp = 0 in X−S 1 −S 2 −S 3 . Consequently, the propagation of the boundary between the infected and uninfected domains in theỸ -direction is suppressed owing to the intended velocities of pedestrians, namely, (±1, 0), as confirmed by comparing of the contour ofp att = 1 in Fig. 5 with that in Fig. 21 . Subsequently, four characteristic infected domains, S 1 , S 2 , S 4 and S 5 , which are defined as S 1 and S 2 are obtained as a result of time evolutions of S 1 and S 2 att = 1, respectively. Consequently, most of pedestrians in S 1 and S 2 recover from the infection, as confirmed bȳ p 1 in S 1 and S 2 . Meanwhile,p in S 4 approximates to its maximum value, namely, unity, att = 1.44. At last, we refer to the local epidemic spread, when intended velocities areṽ 0 i = ±ṽ 0 = ±5, 0, 0 . The number of pedestrians with the intended velocity v 0 , 0, 0 is same as that of pedestrians with −v 0 , 0, 0 . The initial condition is quite same as that in Test A-1. Figure 22 shows contours ofp (upper-left frame),ρ (upper-middle frame) andθ (upper-right frame), whent = 0.168, and contours ofp (lower-left frame),ρ (lower-middle frame) and θ (lower-right frame), whent = 0.96. Certainly, the propagation speed of two characteristic infected domains, which correspond to S 1 and S 2 in the case ofṽ 0 = 1, increases in accordance with the increase inṽ 0 , as shown in upper-left frame of Fig. 22 . Surely, two characteristic infected domains, which correspond to S 4 and S 4 in the case ofṽ 0 = 1, are obtained att = 0.96. Such an increase in the propagation speed inỸ -direction in accordance with the increase inṽ 0 is explained by the fact that the increase inṽ 0 yields the increase in the local temperature. For example,θ ∈ [0.055, 0.095] is obtained, whenṽ 0 = 1 and t = 4.0, whereasθ ∈ [0.83, 3.43] is obtained, whenṽ 0 = 5 andt = 0.96. The increase in the propagation speed of the boundary between the infected and uninfected domains in accordance with the increase inθ is surely confirmed in Fig. 13 . Finally, we confirm that small patterns in contours ofρ andθ appear att = 0.168, whereas such small patterns in contours ofρ andθ change to large patterns att = 0.96. The large patterns in contours of ρ att = 0.96 are self-organized in accordance with the local dense and sparse populations of pedestrians. The large spatial-patterns ofρ andθ att = 0.96 change, temporally, ast increases fromt = 0.96, whereas scales of the large patterns att = 0.96 are not enlarged, ast increases fromt = 0.96. We, however, find that pattern formations in contours ofρ and θ never yield patterns in contours ofp att = 0.96.

A contribution to the mathematical theory of epidemics

Simultaneous reconstruction of evolutionary history and epidemiological dynamics from viral sequences with the birthdeath SIR model

Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models

The mathematics of infectious diseases

Epidemic models and percolation

Velocity and hierarchical spread of epidemic outbreaks in scale-free networks

Immunization strategies for epidemic processes in timevarying contact networks

The role of the airline transportation network in the prediction and predictability of global epidemics

Global disease spread: statistics and estimation of arrival times

Individual-based lattice model for spatial spread of epidemics

Social force model for pedestrian dynamics

Modeling Complex Living Systems: A Kinetic Theory and Stochastic Game Approach

Coarsely grained stochastic Boltzmann equation and its moment equations

Kinetics of influenza A virus infection in humans

Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Claredon

Boltzmann equation for microscopic density

Front propagation into unstable states

A micro-simulation model for pedestrian flows

Hydrodynamics for a granular binary mixture at low density

Front acceleration by dynamic selection in Fisher population waves

Striated populations in disordered environments with advection

A note on the Lattice Boltzmann versus finite difference methods for the numerical solution of the Fisher's equation

Differential Equations and Their Applications

Analytical and numerical study on the nonequilibrium relaxation by the simplified FokkerPlanck equation

Numerical methods for the stochastic Landau-Lifshitz Navier-Stokes equations

Analysis of numerical techniques of dirct simulation Monte Carlo method in the rarefied gas dynamics

Variable soft sphere molecular model for inversepower-law or Lennard-Jones potential

Enskog equations for traffic flow evaluated up to Navier-Stokes order

A fluid-dynamic model for the movement of pedestrians