key: cord-0101337-hvnyusa3
authors: Ha, Seung-Yeal; Park, Hansol; Yang, Seoyeon
title: Relaxation dynamics of SIR-flocks with random epidemic states
date: 2021-10-29
journal: nan
DOI: nan
sha: 8e36a6277be26baadce7e25b8973501a63b31665
doc_id: 101337
cord_uid: hvnyusa3

We study the collective dynamics of a multi-particle system with three epidemic states as an internal state. For the collective modeling of active particle system, we adopt modeling spirits from the swarmalator model and the SIR epidemic model for the temporal evolution of particles' position and internal states. Under suitable assumptions on system parameters and non-collision property of initial spatial configuration, we show that the proposed model does not admit finite-time collisions so that the standard Cauchy-Lipschitz theory can be applied for the global well-posedness. For the relaxation dynamics, we provide several sufficient frameworks leading to the relaxation dynamics of the proposed model. The proposed sufficient frameworks are formulated in terms of system parameters and initial configuration. Under such sufficient frameworks, we show that the state configuration relaxes to the fixed constant configuration via the exponentially perturbed gradient system and explicit dynamics of the SIR model. We present explicit lower and upper bounds for the minimal and maximal relative distances.

The purpose of this paper is to continue the studies begun in [19, 20] on the collective dynamcis modeling of active particles with internal states. Collective behaviors of complex systems are ubiquitous in nature, e.g., crowd dynamics [2] , aggregation of bacteria [8, 9, 10, 26] , flocking of birds [5, 6, 12, 13, 15, 22, 23, 24, 29, 30, 34, 37] , synchronization of fireflies [11, 38] and swarming of fish [16] , etc. See survey articles [1, 3, 7, 17, 21, 33, 36] . In recent years, thanks to the emerging applications to the decentralized control of multi-particle systems, collective behaviors have received a lot of attention from diverse scientific disciplines such as applied mathematics, biology, control theory and statistical physics etc. In this work, we are interested in the first-order modeling of active particles with random epidemic states (susceptible(S), infected(I) and recovered(R)) as an internal state, i.e., we assume that particle i can take aforementioned three epidemic states with certain probabilities denoted by S i , I i and R i whose precise meaning will be clarified in a minute. Thus, the state of particle i is represented by position vector x i ∈ R d and probability vector for epidemic states W i = (S i , I i , R i ) ∈ [0, 1] 3 , respectively. In what follows, we briefly describe how to model the dynamics of position vectors and epidemic vectors via continuous dynamical systems. Let N be the size of system, i.e., the total number of particles in a given ensemble

First, we use the position dynamics x i = x i (t) of the swarmalator model [31, 32] which describes the attractive and repulsive forces:

where κ 2 , κ 3 are nonnegative coupling strengths, and Ψ ij a , Ψ ij r represent attractive and repulsive weights whose explicit functional forms will be discussed in Section 3, and we assume that positive system parameters α, β satisfy the relation:

so that repulsive force is dominant in a small relative distance regime. Here · denotes the standard 2 -norm in R d . Next, we use the modeling spirit of the SIR model for the dynamics of epidemic state W i . We introduce convex set S consisting of all admissible state vectors: and we set S i : the probability that the i-th particle is in susceptible state, I i : the probability that the i-th particle is in infected state, R i : the probability that the i-th particle is in recovered state, and W i := (S i , I i , R i ) : the epidemic probability vector of the i-th particle. Then, we assume that the dynamics of W i is governed by the following coupled system:

Finally, we couple two systems (1.1) and (1.4) via nonnegative system functions Ψ ij a , Ψ ij r , a ij and b i by imposing suitable functional dependences between position variable {x i } and internal variable {W i }:

(1.5) Ψ ij a := Ψ a (W i , W j ), Ψ ij r := Ψ r (W i , W j ), a ij := a( x i − x j ), i, j ∈ [N ]. Moreover, we also assume that there exist positive constants ε a , ε r , M a and M r such that 0 < ε a ≤ Ψ ij a ≤ M a , 0 < ε r ≤ Ψ ij r ≤ M r , i, j ∈ [N ]. See Section 3 for explicit functional relations in (1.5) .

Finally, we combine all three ingredients (1.1), (1.4) and (1.5) together to write down the dynamical system for (W i , x i ) with suitable initial data:

Throughout the paper, we call the coupled system (1.6) as the "SIR-flock model" for simplicity. The goal of this work is to provide sufficient frameworks leading to the emergent dynamics of the SIR-flock model (1.6) . Now, we briefly discuss four main results.

First, we show that if spatial configuration is noncollisional initially, there will be no finite-time collisions so that system (1.6) is globally well-posed by the Cauchy-Lipschitz theory (see Theorem 2.2) .

Second, we present a sufficient conditions for the relaxation of epidemic state W i toward a constant epidemic state. More specifically, our sufficient framework for the relaxation of W i is expressed in terms of network topology (a ij ) and a recovering vector (b i ):

where L is a positive constant and γ is a nonnegative constant. Under this framework (1.7), we show that the epidemic state W i relaxes to a constant state (see Theorem 4.1): there exist a constant state

Third, we show that if initial configuration and system parameters satisfy suitable conditions, then minimal and maximal relative distances are positive uniformly in time, i.e., there exist positive constants δ 1 and δ ∞ such that

Fourth, we show that the spatial configuration relaxes to a constant configuration asymptotically under a suitable condition on the initial configuration and system parameters (see Theorem 6.1).

The rest of this paper is organized as follows. In Section 2, we briefly review basic properties of the SIR epidemic model and the swarmalator model on the asymptotic relaxation of state variables and present a global well-posedness by verifying nonexistence of finite-time collisions. In Section 3, we discuss the modeling spirit of system parameters and coupling functions. In Section 4, we study the relaxation of epidemic states toward a constant state and in particular, we provide a sufficient condition leading to asymptotic removal of infected particles. In Section 5, we study the existence of positive lower bound and upper bound for the minimal and maximal relative distances, respectively. In Section 6, we study a relaxation of spatial configuration toward a fixed spatial configuration using the perturbed gradient flow theory. In Section 7, we provide several numerical examples and compare them with analytical results obtained in previous sections. Finally, Section 8 is devoted to a brief summary of our main results and some remaining issues for a future work.

In this section, we briefly review basic properties on two related models "the SIR epidemic model" and "the swarmalator model" which correspond to the subsystems of system (1.6).

2.1. The SIR epidemic model. In this subsection, we briefly discuss the SIR model [4, 35] which is a prototype model for the spread of disease or virus. Kermack and McKendrick's work in 1927 motivated a large number of modelings in epidemics [25] . In particular, there are some of the works for diffusive SIR models on stability and numerical methods [18, 14, 28] . First, we introduce the following three observables: S = S(t) : ratio of susceptible individuals in a total population, I = I(t) : ratio of infected individuals in a total population, R = R(t) : ratio of recovered individuals in a total population.

Then, we assume that the state (S, I, R) is governed by the Cauchy problem to the coupled system of ODEs:

where a and b are positive constants. This model was designed to explain the spreading of epidemic diseses. We can express the SIR model via the following pictorial diagram.

Next, we list basic properties of the SIR model in the following lemma.

Proposition 2.1. Let (S, I, R) be a global solution of (2.1) with initial data satisfying

Then, one has

Proof. (i) We add three equations in (2.1) to geṫ S +İ +Ṙ = 0.

This yields

(ii) We integrate the first two equations in (2.1) to obtain

3)

The nonnegativity of R follows from the third equation and the nonnegativity of (2.3) 2 :

On the other hand, it follows from (2.2) and (2. 3) that R is bounded by 1: 

2.2. The swarmalator model. Let x i and θ i be the position and phase of the i-th swarmalator, respectively. Then, its dynamics is governed by the Cauchy problem to the Kuramoto-type swarmalator model:

Here, constants w i and ν i are the natural velocity and frequency of the i-th particle, respectively, · denotes the standard Euclidean 2 -norm in R d . System parameters and interaction functionsΨ a andΨ r satisfy

If we simply set γ = 0, system (2.5) 2 becomes the Kuramoto model [27] :

Due to the singular terms in the swarmalator model (2.5) 1 , it is important to make sure that there is no finite-time collision which can be quantified in the following proposition.

First, we set several parameters:

(2.7)

Next, we state two results on the positivity of minimal and maximal distance between particles.

Proposition 2.2. (Positivity of minimal distance) [20] Suppose system parameters satisfy (2.6), and the initial data (X 0 , Θ 0 ) is non-collisional:

Then, there exists a global solution (X, Θ) to (2.5) -(2.6) andδ 1 > 0 such that

[20] Suppose system parameters and initial data satisfy the following framework:

• (F1) : For α = 1,

where y * is the largest root of the equation:

and let (X, Θ) be a solution to (2.5) -(2.6). Then, there exists a positive constantδ ∞ such that sup 0≤t<∞ D(t) ≤δ ∞ .

Proof. For a proof, we refer to [20] .

Consider the following ansatz forΨ a andΨ r : Ψ a (θ) := 1 + J cos θ,Ψ r := 1, |J| < 1.

In this case, system (2.5) becomes (2.10)

By Proposition 2.1 and convergence result of the perturbed gradient system, system (2.10) exhibits phase synchronization.

Theorem 2.1.

[20] Suppose system parameters, natural velocities and initial data satisfy one of the frameworks (2.8) and (2.9), and κ > 0, D(ν) = 0 and 0 < D(Θ 0 ) < π, and let (X, Θ) be a solution to (2.10) . Then there exists asymptotic state X ∞ such that

Proof. For a proof, we refer to Section 4.2 of [20] .

A global well-posedness. In this subsection, we discuss a global well-posedness of system (1.6). Since the R.H.S. of (1.6) 2 contains the term x j − x i in the denominators, as long as we can rule out the possibility of finite-time collisions, we obtain a global wellposedness using the standard Cauchy-Lipschitz theory. In the sequel, we show that system (1.6) does not admit a finite-time collision, as long as there is no collisions initially. This will be done using a contradiction argument and Gronwall's inequality.

Suppose that initial spatial configuration satisfy min 1≤i =j≤N

Then, by the continuity of solution, there will be no collisions between particles at least small time interval [0, ε) with ε 1. Then, using the Cauchy-Lipschitz theory, we can show that system (1.6) has a local smooth solution {(W i , x i )} in the time-interval [0, ε). Now, in order to show a global well-posedness, it suffices to show that there will be no finite-time collisions. Suppose there is a finite-time collision and let t 0 be the first collision time. We take one of the particles that make a collision and fix it by p. Then, we define a set containing all the particles involved in the collision at time t 0 by C. Now, we set the following handy notation:

Then, it is easy to see that

By direct calculations, one has

(2.12)

In the following lemma, we provide estimates for I 1i .

where |A| is the cardinality of the set A and δ := inf i∈C, j∈R

Proof. We provide estimate for I 1i (i = 1, 2, 3, 4) as follows:

(i) (Estimate of I 11 + I 13 ): Note that

Therefore, one has (2.13)

Similarly, one has (2.14)

We combine (2.13) and (2.14) to obtain

(ii) (Estimate of I 12 + I 14 ): Since |x k − x i | ≥ δ, one has

where we used the Cauchy-Schwarz inequality to see

Similarly, one has

It follows from (2.12) and Lemma 2.1 that

where c 1 , c 2 , and c 3 are positive constants defined as follows.

On the other hand, since 1 ≤ α < β in (1.2), there exists a small positive constant ε 1 such that for t ∈ (t 0 − ε, t 0 ),

Therefore, we combine (2.15) and (2.16) to find

Moreover, there exists a positive constant ε * ∈ (0, ε) such that for t ∈ (t 0 − ε * , t 0 ),

. We integrate the above differential inequality from t 0 − ε * to t 0 to get

. On the other hand, it follows from (2.11) that the left-hand side is zero, but the right-hand side is strictly positive, which is contradictory. Therefore, there are no finite-time collisions between particles from the well-prepared initial configuration. Finally, we can summarize the previous argument as follows.

Theorem 2.2. (A global well-posedness) Suppose that initial spatial configuration is noncollisional in the sense that min

1≤i =j≤N

Then, there exists a unique global solution {(W i , x i )} to system (1.6) in any finite-time interval.

In this section, we discuss meaning of the system parameters and system functions appearing in the SIR-flock (1.6):

(a ij ) : interaction topology, (b i ) : recovering vector, Ψ a , Ψ r : coupling weight functions.

3.1. Interaction topology and recovering vector. For the modeling purpose, we assume that the nonnegative value a ij = a ij ( x i (t) − x j (t) ) depends on the relative distance x i − x j between the i-th particle and the j-th particle. To be definiteness, we set

where L is a positive constant and γ is a nonnegative constant. Since the i-th particle can not be affected by itself, we assume the diagonal entries of (a ij ) are zero:

Now we discuss the natural recovering vector b = (b 1 , · · · , b N ). Suppose that every particles have own immune system, the disease will disappear automatically. Thus, we

. If the immune system of the i-th particle is well-functioning, then b i has a large value. In contrast, if the immune system of the i-th particle is not well-functioning, then b i mpimwwwi will have a small value. In this work, we assume that the natural recovering vector b i is a positive constant (see the following figure in the sequel).

Coupling weight functions. Recall that the dynamics of x i is governed by the following system:

Now, the matter of question is how to model Ψ a and Ψ r in terms of W i s.

If two state vectors W i and W j are similar, then particles x i and x j will attract each others, whereas if two state vectors W i and W j are dissimilar, then x i and x j will repel each other. For definiteness, we set

where α and β are positive constants. Here ε a and ε r are positive constants presenting social distancing. These social distancing constants play an important role in preventing collisions between particles. Note that the ansatz (3.3) satisfies symmetry:

Next, we will impose conditions on α and β later in (3.2) . To illustrate the functional relations in the right-hand side of (3.3), we consider an ensemble which is partitioned into two sub-ensembles:

{Infected particles} or {susceptible or recovered particles}.

Recall that I i and S i + R i are probabilities that i-th particle are in infected state and are not in infected state, respectively.

Consider the inner product-like function G as follows:

Then, Ψ a and Ψ r satisfy the following monotonicity properties:

(1) If W i and W j become similar, Ψ a increases and Ψ r decreases.

(2) If W i and W j become dissimilar, Ψ a decreases and Ψ r increases.

The relations (3.5) and (3.7) yield (3.3). On the other hand, it follows from (3.3) and (3.6) that the coupling weight functions Ψ a and Ψ r admit positive lower bound and upper bounds:

Proposition 3.1. Suppose that system parameters satisfy

and let (S, I, R) and (S i , I i , R i ) be the solutions of system (2.1) and system (1.6) with the initial data:

Then we have

Proof. Suppose the conditions (3.9) hold. Then, the relation (3.1) and (2.1) become

Since (S i , I i , R i ) and (S, I, R) satisfy the same system (2.1) with the same initial data, by the uniqueness of ODE solution, we have the desired uniqueness (3.10).

3.3. Symptom expression vector. Next, we introduce the symptom expression vector s = (s 1 , · · · , s N ) motivated by the ongoing pandemic COVID-19. In April 2020, the daily COVID-19 confirmed number per day of Korea was less than 20. Most of confirmed people were isolated, however patients with no symptoms of COVID-19 were still spreading the virus. That is why the daily confirmed number per day does not goes to 0 directly. So we need to consider the ratio of symptom expressions. Some people show symptoms well when they got the virus, in contrast some people does not show any symptoms although they got the virus. So even for the same inputs, different degree of output can emerge. Thus, we define the symptom expression vector that are related to previous phenomena. We define define each component s i ∈ [0, 1] to following the following properties:

• If the i-th particle shows the symptom well, then we put large value to s i .

• If the i-th particle does not show the symptom well, then we put small value to s i .

We may use the symptom expression vector s to express the condition of being suspected as a patient. In this paper, we assume that the symptom expression vector s is a timeindependent constant vector. We assume that the condition of being suspected as a patient only depends on the product s i I i , and we also define some fixed threshold constant c to make decision. If s i I i (t) ≥ c, then the i th person will be confirmed at time t.

In this section, we study the relaxation dynamics of the SIR-flock model (1.6). First, we show that S is invariant along (1.6). 

where S is the set defined in (1.4) .

Proof. We basically use the same arguments as in Proposition 2.1.

(i) (Verification of the second relation in (4.1)): It follows from (1.6) 1 that

This yields

Now, we need to show the positivity of S i , I i and R i .

• (Positivity of S i ): Due to (4.2), it suffices to check the positivity of S i , I i and R i . Now, we integrate (1.6) 1 to find

• (Positivity of I i ): It follows from (2.1) 2 thaṫ

This yields Since each I i (t) is analytic, there exists the refinement of time-interval 0 = t 0 < t 1 < t 2 < · · · such that the index of I mt is not changed on each interval [t j , t j+1 ), and

Therefore, we have

This yields

Therefore, one has

• (Positivity of R i ): Since S i + I i + R i = 1, S i ≥ 0 and I i ≥ 0, one has

Therefore, we combine all the estimates for S i , I i and R i to get

(ii) (Verification of the second relation in (4.1)): Recall the equation for x i :

Since the right-hand side of the above equation is skew-symmetry with respect to interchange map i ←→ j using (3.4), the total sum i x i is time-invariant. Now, we discuss the asymptotic convergence of the probability vector W i toward a constant probability vector W ∞ ∈ S in (1.4). (1) There exists a constant state

(2) If the recovering value b i satisfies

then there exists a positive constant λ := min

Proof. (i) Since I j and S i are non-negative,

which means S i is non-increasing. Since S i is bounded below by zero, there are some constants S ∞ i that S i (t) converges to S ∞ i , as t goes ∞ by monotone convergence theorem. Similarly,Ṙ i = b i I i ≥ 0, and R i is non-decreasing and bounded above, there are some constants R ∞ i that R i (t) converges to R ∞ i , as t goes ∞. Moreover, as R i converges, 0 = lim (ii) Recall the equation for I i :İ

This implies

Next, under a more relaxed condition compared to (4.3), we improve the second result of Theorem 4.1 as follows. 

then N i=1 I i decays to zero exponentially fast. Proof. First, we write the dynamics (1.6) 1 of I i s in matrix form:

where I and A are given as follows.

I(t) := (I 1 (t), · · · , I N (t)) and A(t) := (a ij (t)) 1≤i,j≤N .

Since each of S i is decreasing,

This implies that for every t,

where ≤ is a partial order compoentwise, and

Therfore, we obtain (4.5)İ = κ 1 diag(S 1 , S 2 , · · · , S N )AI − bI ≤ (κ 1Ā − bI)I .

We multiply the both sides of (4.5) by 1 to yield

where λ i is given as follows.

the condition (4.4) gives a more relaxed condition on {b i } compared to (4.3).

In this section, we provide explicit quantitative estimates on relative distances between particles. For a given spatial configuration {x i }, we consider the set of all relative distances { x i − x j } for i = j, and rearrange them in an increasing order:

In what follows, we derive a uniform lower-bound for D 1 and a uniform upper-bound for D Q so that we can control the singular terms in (1.6) uniformly in time. First, we recall notation: (1) (Existence of a positive lower bound to D 1 ): If initial data satisfy min 1≤i =j≤N

there is a positive constant δ 1 > 0 such that

(2) (Existence of a positive upper bound to D Q ): If initial data and system functions satisfy

there exists a positive constant δ ∞ such that

Proof. We leave its proof in the next two subsections.

A positive lower bound for minimal relative distance. In this subsection, we show that D 1 (t) has a positive lower bound δ 1 which is independent of t. Since the particles do not collide in finite-time interval, there exists an analytic solution in finite time (see Theorem 2.2). Therefore, each distance D ij and ordered distance D i are Lipschitz continuous. In addition, according to the analyticity, for given D ij , there exist a sequence of times (t n ):

such that we can decompose the whole time interval [0, ∞) into the union of subintervals

to make sure that D ij is not changed in each subinterval. To demonstrate the existence of positive uniform lower bound, we first provide three lemmas. First, we show a positive lower bound of maximal distance D Q (t):

Next, we provide a series of lemmas to estimate maximal and minimal relative distances. First, we derive a differential inequality for D.

Proof. For given l ∈ N ∪ {0}, we choose i and j such that

Then, by direct calculation, we obtain

for t ∈ T l . We now estimate I 2i (i = 1, 2, 3) one by one.

• (Estimate of I 21 ): We use

• (Estimate of I 22 ) : Since x i − x j is the maximal distance,

Similarly, we have

Therefore, one has

(5.5)

• (Estimate of I 23 ): Similar to the estimate of I 22 , we obtain To derive a positive lower bound for D 1 , we first verify that the maximal distance D is bounded away from zero uniformly in time, and then we show that this positive lower bound propagates to the lower graded relative distance when we reach to the minimal relative distance D 1 (the first assertion in Theorem 5.1) via the following two lemmas. Now we derive a positive lower bound for the maximal distance D(t).

Lemma 5.2. (maximal relative distance) Suppose that initial spatial configuration is noncollisional: min

1≤i =j≤N

Proof. We use an induction argument on the time intervals T l = [t l−1 , t l ).

• (Initial step): we claim that

Suppose not, since δ Q ≤ D(0), there exist t 11 and t 12 such that 0 ≤ t 11 < t 12 < t 1 , D(t 11 ) = δ Q and D(t) < δ Q in (t 11 , t 12 ]. Therefore, for t ∈ (t 11 , t 12 ],

This implies

which is contradictory to D(t) < δ Q in (t 11 , t 12 ].

• (Inductive step) : Suppose that for some l ≥ 1,

Now we consider D(t) on T l+1 . Due to the continuity of D(t) and induction hypothesis, one has

Therefore, we can use the same criteria above to get

Thus, one has the desired estimate.

In the following lemma, we show the backward propagation of a positive lower bound for the ordered relative distance. Proof. We refer to the proof of Lemma 3.4 in [20] . Now, we are ready to provide a positive lower bound for the minimal relative distance D 1 .

Proof of the first assertion in Theorem 5.1: Suppose initial spatial configuration is noncollisional:

and let {(W i , x i )} be a solution to system (1.6). By Lemma 5.2, there exists δ Q > 0 such that

For q = Q, we apply Lemma 5.3 to show that

Again, we apply Lemma 5.3 inductively until we reach q = 1 to derive the desired a positive lower bound for D 1 .

A positive upper bound for maximal relative distance. In this subsection, we derive an existence of a positive upper bound for the maximal relative distance. First, we study the estimate for a differential inequality.

Lemma 5.4. Let y : [0, ∞) → R + be a differentiable function satisfying the following differential inequality:

where constants a, b, p and q satisfy

Then, there exists a positive constant δ such that sup 0≤t<∞ y(t) ≤ δ.

Proof. We use phase line analysis forż = −az −p + bz −q . For this, we define

Then, it is easy to see that f (y) > 0 for y < y * , f (y * ) = 0 and f (y) < 0 for y > y * . Therefore, y is uniformly bounded by max{y 0 , y * }. Now, we are ready to provide a proof for the second assertion in Theorem 5.1.

Proof of the second assertion in Theorem 5.1: Suppose initial data and system parameters satisfy

and let {(W i , x i )} be a solution to system (1.6). Then, for t ∈ (0, ∞), we choose i = i t and j = j t such that

In what follows, we estimate I 3i (i = 1, 2, 3) one by one.

• (Estimate of I 31 ): We use

• (Estimate of I 32 ): We use D(t) = x i − x j to see

Now, we define the angle between x i − x j and x k − x i by θ ijk . Since

one has cos(θ ijk ) ≤ 0. Thus, one has

Since

we have (5.9) I 32 ≤ 0.

• (Estimate of I 33 ): Similar to I 32 , we have (5.10)

Finally, in (5.7), we combine all the estimates (5.8), (5.9) and (5.10) to find d dt

to derive the desired estimate.

Remark 5.1. Note that for 1 ≤ α < β, nonexistence of finite-time collisions and uniform lower bound of diameter are always warranted, however conditions on parameter is needed to guarantee uniform upper bound of diameter.

In this section, we study the relaxation of spatial configuration toward a constant configuration.

Consider a perturbed gradient system with decaying forcing f :

where V is a one-body potential. Proposition 6.1. Suppose that external forcing f and potential V satisfy the following conditions:

(1) There exist positive constants C and λ such that

There exists a compact set K such that x(t) is contained in K for all t ≥ 0, (3) V is analytic in K and the map t → |∇ x V (x(t))| 2 is uniformly continuous in time, and let x = x(t) be a solution of (6.1). Then, there exists x ∞ ∈ K such that

Proof. We refer to [20] for a proof. Now we are ready to prove the convergence of spatial configuration. Theorem 6.1. (Relaxation of spatial configuration) Suppose initial data and system parameters satisfy

where Λ is a positive constant defined in (5.1), and let {(W i , x i )} be a solution to system

Proof. We consider three cases below.

• Case A (α = 2 and β = 2): We set the potential function V and the function f bẏ

Note that X is in R N d when each of x i is d−dimensional. First of all, we claim that the analyticity of the potential function V defined in R N d . we can restrict the solution space because of the existence of uniform lower bound of distance. We may exclude the hyperplanes E ij := {x ∈ R N d : 

By the continuity of the solution, O k are invariant sets for X. Without loss of generality, we may assume X(0) ∈ O 1 and thus X(t) ∈ O 1 . That is O 1 is an open domain and V is analytic on O 1 . Now, we will show that there exists positive constants C and λ such that

Note that we have verified that there exists a positive constant c such that

It's clear that for all i, I i (t) ≤ N e −ct =: w.

In addition, one has

This means Ψ a (W i , W j ) − 1 − ε a and Ψ r (W i , W j ) − ε r decay in exponential order. Since we have uniform upper and lower bounds of x i − x j , there exist positive constants C and λ such that

Lastly, upper and lower bounds guarantee the boundedness of V X (X(t)) and ∇ X V (X(t)). Similarly, d dt |∇ X V (X(t))| 2 is bounded too. Thus we get uniform continuity of |∇ X V (X(t))| 2 in time. Thus, every hypothesis in Proposition 6.1 is satisfied to get X ∞ that X(t) converges.

• Case B (α = 2): By setting

we can prove the desired estimate similarly as in Case A.

• Case C (β = 2): Similar to Case B, we can show the desired estimate with

Next, we consider a two-particle system with the following initial SIR state:

S 1 (0) = s 1 , I 1 (0) = 1 − s 1 , R 1 (0) = 0, S 2 (0) = s 2 , I 2 (0) = 1 − s 2 , R 2 (0) = 0, and we set s := max{s 1 , s 2 }.

Consider a system for spatial position:

Note that sinceẋ 1 +ẋ 2 = 0, the center of the two particles is constant:

and derive the equation for x:

As long as each coordinate of x is positive,

Thus, one has

so the difference x has upper and lower bounds that are uniform in time.

It follows from Lemma 5.2 that for δ Q = min D(0), κ 3 mr

and D(t) = x(t). Therefore, x(t) has both upper and lower bounds and by the same reason in Theorem 6.1, there exists x ∞ such that x(t) converges, as t tends to infinity.

This yields Ψ 12 a = 1 + ε a and Ψ 12 r = ε r . Sinceẋ ∞ = 0, we have

In addition, we can weaken the condition that guarantees the exponential decay of I i .

Theorem 6.2. Suppose system parameters satisfy

and let {(W i , x i )} be a solution to system (2.1). Then, there exist positive constants A, ε * and t 0 such that

Thus, for every ε > 0, there exists positive constant t ε with respect to ε such that

By the given condition, there exists a positive constant ε * such that

For t > t ε * ,

This yields

Remarks. Since

The Jacobian matrix of system at (I 1 , I 2 ) = (0, 0) is

Next, we consider two cases for b.

• Case A: Consider the case

In this case, all eigenvalues are real and eigenvalues for I ! and I 2 are negative. This means that I 1 and I 2 decay to (0, 0) exponentially fast.

The right-hand side is less or equal to

Note that the left-hand side ∞ 0 (I 1 + I 2 )(t)dt denotes the number of total infections. Thus, we can say that I i (t) decays to zero exponentially fast, or the total number of total infection is bounded.

In this section, we provide several numerical examples to confirm analytical convergence results that we have shown in previous sections. At first, we set initial condition by (S i (0), I i (0), R i (0)) = (1, 0, 0) for uninfected 16 people and (S i (0), I i (0), R i (0)) = (0.1, 0.9, 0) for infected 4 people. Initial location is set by random seed in 3 by 3 plane. We set κ 1 = 1, ε a = ε r = 0.2, and we use the fourth order Runge-Kutta method for all simulations.

• (Convergence of (S i , I i , R i )): In Figure 1 and Figure 2 , we used b = 0.4, γ = 1, L = 1, here. In Theorem 4.1, we have shown that S i , I i , and R i converges. To authenticate this assertion by numerical method, we observed convergence of S i , I i , and R i . Moreover, by (a) S, I, R for α = 1, β = 2, κ2 = 1, κ3 = 5 (b) S, I, R for α = 2, β = 5, κ2 = 10, κ3 = 1 Figure 1 . Convergence of S, I and R observingṠ i ,İ i , andṘ i that converges to zero, we obtain that S i , I i and R i converge. and observed two simulations based on two set of system parameters:

(α, β, κ 2 , κ 3 ) = (1, 2, 1, 5), (2, 5, 10, 1).

(a) S, I, R for α = 1, β = 2, κ2 = 1, κ3 = 5 (b) log i Ii for α = 1, β = 2, κ2 = 1, κ3 = 5 Figure 3 . Exponential decay of log i I i • (Effect of social distancing): Since κ 2 and κ 3 are coefficients of attracting and repulsion forces, respectively, larger ratio κ 3 κ 2 represents for more intensive social distancing. We observe that the maximal value of i I i decreases, as κ 3 κ 2 increases. We set L = 1 and κ 2 = 1, and changed values of b and Ψ with κ 3 = 1, 5, 10, 50. We plot graph for the cases of κ 3 κ 2 = 1, 5, 10 and 50. • (Behavior of the particles): We observed that infectious particles aggregated and they were isolated from non-infectious factors. It implies that this model can effectively reduce the number of infectious particles by adjusting the coefficients. We illustrated the behavior of the particles in b = 0.4, γ = 3, L = 1, and κ 2 = 1, κ 3 = 10. The purple ones represents for non-infectious particles. Moreover, we compared with the model that has no attracting and no repulsing term by setting κ 2 = κ 3 = 0. It is illustrated in the first graph in Figure 6 , red one is for κ 2 = 1, κ 3 = 10 and blue one is for κ 2 = κ 3 = 0. We could decrease the average of I i (t). We set b = 0.2, γ = 3, L = 3, α = 2, β = 5 and repeated the process. In this paper, we have studied the emergent behaviors of a flock with SIR internal states, and have presented a new particle model with an aggregate property and epidemic internal forces by combining the SIR model and the swarmalator model. We considered that each particle has a state vector which is a kind of internal state following the SIR dynamics. From this argument, we could obtain the SIR-aggregation model by imposing the repulsive/attractive force between two particles depending on the distance between two particles and the internal states. We imposed strong repulsive force if one of them has a high probability of infection to model the social distancing, in contrast, we imposed weak repulsive force if both of them has low probabilities of infection. From this modeling, we modeled the social distancing of particles. We expect that we can do numeric experiments to find the efficiency of social distance for given strategies. We also provided the theoretical result of the SIR model, for example, nonexistence of finite-time collisions, positive lower bound for minimal relative distance, and a uniform upper bound for spatial diameter. From the numerical simulations, we could check that the isolation of the infected particle is an effective way to reduce the number of infected particles. Since we only considered that the recovered particles were never infected again, we have monotonicity on the number of susceptible/recovered particles. From those properties, we could prove the emergent dynamics. In our proposed model, we did not consider the reinfection which can happen in reality. Thus, we will leave this issue for a future work. 

The Kuramoto model: A simple paradigm for synchronization phenomena

Towards a systems approach to behavioral social dynamics

Vehicular traffic, crowds, and swarms: From kinetic theory and multiscale methods to applications and research perspectives

Analysis of SIR epidemic model with information spreading of awareness

Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study

A quest toward a mathematical theory of the dynamics of swarms

Towards a mathematical theory of behavioral swarms

Finite-time blow-up of L ∞ -weak solutions of an aggregation equation

Blow-up in multidimensional aggregation equations with mildly singular interaction kernels

L p theory for the multidimensional aggregation equation

Biology of sychronous flashing of fireflies

Sharp conditions to avoid collisions in singular Cucker-Smale interactions

Flocking and turning: a new model for self-organized collective motion

Numerical modelling of an SIR epidemic model with diffusion

Emergent behavior in flocks

Large-scale dynamics of the Persistent Turing Walker model of fish behavior

A kinetic flocking model with diffusion

Localized outbreaks in an S-I-R model with diffusion

A mean-field limit of the particle swarmalator model

Emergent behaviors of the swarmalator model for position-phase aggregation

Collective synchronization of classical and quantum oscillators

A simple proof of Cucker-Smale flocking dynamics and mean-field limit

Emergent dynamics of a thermodynamically consistent particle model

From particle to kinetic and hydrodynamic description of flocking

A contribution to the mathematical theory of epidemics

Emergent behaviour in multi-particle systems with non-local interactions

International symposium on mathematical problems in mathematical physics

Analysis of a local diffusive SIR model with seasonality and nonlocal incidence of infection

Heterophilious dynamics enhances consensus

A new model for self-organized dynamics and its flocking behavior

Ring states in swarmalator systems

Oscillators that sync and swarm

Synchronization: A universal concept in nonlinear sciences

Flocks, herds, and Schools: A quantitative theory of flocking

On the global stability of SIS , SIR and SIRS epidemic models with standard incidence

Collective motion

Novel type of phase transition in a system of self-driven particles

Biological rhythms and the behavior of populations of coupled oscillators