key: cord-0176886-fkya7les
authors: Kobayashi, Hisashi
title: Stochastic Modeling of an Infectious Disease Part II: Simulation Experiments and Verification of the Analysis
date: 2021-01-26
journal: nan
DOI: nan
sha: 8a33b94dbc27beabdd228acf3ec2d7700550b1e9
doc_id: 176886
cord_uid: fkya7les

In Part 1, we introduced a stochastic model of an infectious disease, based on the BDI (birth and death with immigration) process. We showed that random processes defined by this model can capture the essence of the stochastic, often erratic, behavior of the infection process. The most significant finding was that it is negative binomial distributed with small r, hence it is geometrically distributed with an exceedingly long tail. This leads to a much larger disparity in the epidemic patterns than has been known to the modeling community. In Part 2 we conduct simulation experiments by implementing an event-driven simulator. Several independent runs are presented to verify the findings reported in Part 1. The enormous variations among the sample paths obtained from several consecutive and independent runs with statistically identical conditions confirm our analysis. Note that the epidemic pattern that we observe is merely one sample path taken out of infinitely many paths. Thus, one such path cannot represent the ensemble of the underlying random process. By plotting the simulation runs in the semi-log scale, we can see that the probabilistic chances in the early phase of the infection process determine the behavior of the process. Once the process has reached a considerable number, the weak law of large numbers sets in, and the process behaves less erratically than in the early phase. One important implication of our findings is that it would be a futile effort to attempt to identify all plausible causes of the epidemic patterns. Mere luck may play a more significant role than most people believe. It would be worth remarking that the BDI process might be applicable to explain the disparity in wealth among individuals of similar earning power, expenditure, and investment portfolio.

where I 0 is the initial value:

and 1 . A(t) is the cumulative count of infected arrivals from the outside. We assume the arrival pattern is completely random, i.e., a Poisson process with rate ν [persons/day].

2. B(t) is the cumulative count of the infections that occur in the interval [0, t]. An infection occurs at the rate λ [person/day/infectious person].

3. R(t) is the cumulative count of the infected persons who recover, are removed or die in [0, t].

This event occurs at the rate µ [persons/day/infected person].

The P n (t), n = 0, 1, 2, · · · should satisfy the following set of differential equations:

dP n (t) dt = [(n − 1)λ + ν]P n−1 (t) − [n(λ + µ) + ν]P n (t) + (n + 1)µP n+1 (t), n = 1, 2, 3. · · · , (5) with the initial condition (3) .

In order to find solve the above differential equations, we use the probability generating function (PGF) defined by

The set of differential equations (4) and (5) are then transformed into one partial differential equation (PDE):

∂G(z, t) ∂t = (z − 1) (λz − µ) ∂G(z, t) ∂z + νG(z, t) ,

with the condition G(z, 0) = z I 0 .

We apply Lagrange's method to the above PDE, obtaining

where a λ − µ, and r ν λ .

If I 0 = 0, (9) reduces to G(z, t) = a λz − µ − λ(z − 1)e at

Copyright ©Hisashi Kobayashi 2021 where β(t) λ(e at − 1) λe at − µ . (12) To obtain E[I(t)] = I(t), we use the formula E[I(t)] = ∂G(z,t) ∂z z=1

:

It would be easy to observe that

• If a > 0, I(t) grows exponentially to infinity, as t → ∞;

• If a < 0, it converges to − ν a > 0. • If a = 0, I(t) = νt for all t ≥ 0.

Let T [days] be the number of days that is required for I(t) to double. Then from (13) , we find e aT = 2, or aT = ln 2 = 0.693.

If a = 0.2 [day −1 ] (as in the running example of Part I and the present paper), I(t) doubles in every T ≈ 3.5 [days], hence quadruples every week. 1 In order to obtain the PMFs P n (t), by referring to the definition of the PGF (6), we differentiate G(z, t) (11) w.r.t. z, set z = 0, divide the resulting expression by n!. That is, P 0 (t) = G(0, t) = (1 − β(t)) r P 1 (t) = ∂G(z, t) ∂z z=0 = (1 − β(t)) r rβ(t) P 2 (t) = 1 2! ∂ 2 G(z, t) ∂z 2 z=0 = (1 − β(t)) r (r + 1)r 2! β(t) 2

. . . P n (t) = 1 n! ∂ n G(z, t) ∂z n z=0 = (1 − β(t)) r (n + r − 1) · · · (r + 1)r n! β(t) n ,

which we can write as P n (t) = K(n, r)(1 − β(t)) r β(t) n , n = 0, 1, 2, · · · ,

where K(n, r) is the binomial coefficient:

K(n, r) = (n + r − 1)(n + r − 2) · · · (r + 1)r n! n + r − 1 n .

The above distribution P n (t) takes the form of a Pascal distribution, or negative binomial distribution (NBD), often denoted as NB(k, q), which is the probability of the number of failures n in Bernoulli trials needed to achieve k successes, where q is the probability of failure per trial P N BD n = K(n, k)(1 − q) k q n , n = 0, 1, 2, · · · .

The PGF of this NBD is given by

Thus, by equating k = r ν λ and q = β(t), we see from (11) that the BDI process I(t) is negative binomial distributed according to NB(r, β(t)) at given t.

We note the following properties of K(n, r) (17) and the NBD (16):

i K(0, r) = 1, thus P 0 (t) = (1 − β(t)) r for all t ≥ 0.

ii K(1, r) = r, thus P 1 (t) = r(1 − β(t)) r β(t). Therefore, P 1 (t) = rβ(t)P 0 (t).

iii K(n, 1) = 1 for all n. Then the PMF reduces to a geometric distribution (1 − q)q n with q = β(t)

iv K(n, 2) = n + 1 for all n;

v If r < 1, then K(n, r) < 1 for all n ≥ 0. From Property ii and (16), we find that for each t, the PMF P n (t) is a monotone decreasing function of n, i.e., P n (t) ∝ β(t) n , where β(t) approaches 1 from below as t becomes large.

vi If 1 < r < 2, then K(n, r) is a monotone increasing concave (i.e., convex cap) function of n, bounded by 1 < K(n, r) < n + 1.

vii If r > 2, then K(n, r) is a monotone increasing convex function of n, bounded from below K(n, r) > n + 1.

viii By using Stirling's approximation for factorials and the Gamma function,

we can approximate K(n, r), for n > 1 and r > 1, by

Thus, if r is small, e.g., r ≈ 1, the distribution of (16), for any given t, slowly decreases towards zero as n increases. Thus, the distribution has a long tail (see Figures 7-12 of Part I). Consequently, different sample paths of I(t) are expected to exhibit enormous disparities, as will be shown in the next section, where we present simulation results, which demonstrate huge spreads across different sample paths.

It goes without saying that the most decisive factor that shapes the epidemic pattern is the exponential parameter a = λ − µ. However, for a given a, different simulation runs show enormous differences in their infectious patterns. This wild random behavior is perhaps well beyond what most epidemiologists are cognizant of. In other words, we should accept that we won't be able to find out all factors that can explain why some cities or countries are experiencing more hardship than others in terms of the number of infections and death tolls.

This also implies that Japan and other countries that have observed a relatively small number of infections and casualties should be aware that they may not be as fortunate in new waves of the pandemic. They should be well prepared for the worst case scenario in order to protect their citizens.

In order to be able to provide specifically what will be the worst possible scenario, we need to come up with an accurate estimates of the model parameters λ, ν and ν, by carefully analyzing real and reliable data, and make the most likely estimates (or ranges of estimate), which will be the main focus in Part IV [5] of our report.

Simulations are often used when analytic techniques to estimate or predict the performance of a complex system are hard to come by. So-called Monte-Carlo simulation techniques, such as variance reduction techniques, are often adopted to estimate a given performance measure accurately and efficiently. The main purpose of our simulation experiments here, however, is different from these situations in that we fortunately have obtained an exact analytic result of the system, i.e., we have found closed-form expressions for the PMFs P n (t)'s by solving a partial differential equation.

We present here the results of various simulation runs primarily to better explain and demonstrate the validity and utility of our stochastic model. Plots of various sample paths of our simulation experiments should serve as a convincing and understandable evidence to support our model.

We shall briefly describe how our simulator is designed for the benefit of some readers who may not be sufficiently familiar with probabilistic simulation 2 of a Markov process such as the birthdeath-immigration process. We adopt the time-asynchronous approach as opposed to the timesynchronous approach, which may be more suited to simulating events that occur periodically at predetermined moments, such as seen in time-synchronous communication systems. The timeasynchronous approach is also referred to as the event-scheduling approach.

In our model, there are three types of events; (A) an arrival of an infected person from the outside; (B) an infection caused by an infectious person; and (R) a recovery/removal/death of an infected person.

The random process I(t) defined by the differential equations (4) and (5) is a Markov process. At an arbitrarily given time t, when I(t) = n, the time interval until the next event is a random variable (RV) X with a negative exponential distribution of mean 1/γ(n), 3 where

That is, the distribution function of the RV X is

Generation of instances of the random variable X can be done by calling the random number generator (RNG) that is available in almost all programming languages, including MATLAB. 4 It will generate an instances u of the random variable U , which is uniformly distributed between 0 and 1. Then we transform u by the function F −1 X (·), i.e., find x such that u = F X (x) = 1 − e −γ(n)x . In other words, x = log(1 − u)/γ(n). Since u is uniformly distributed between 0 and 1, so is 1 − u, thus x log u/γ(n) can be used instead of x, At the occurrence of an event, we classify it as one of the three types of events as follows: choose type-(A) with probability P A ν/γ(n); choose type-(B) with probability P B = nλ/γ(n); and choose type-(R) with probability P R = 1 − P A − P B = nµ/γ(n). Note that if n = 0, then P A = 1, and P B = P R − 0.

The rationale for this extremely simple event scheduling approach is that the Poisson process possesses the following beautiful properties: (a) the memoryless property, (b) the reproductive additivity, and (c) the decomposition property. 5 Thus, the core of our simulation program is as follows.

1. If I(t) = 0, then γ(0) = ν, hence, the only possible event to consider is a type-(A) event.

The time until the next arrival, x, is determined by x = log u/ν, where u ∈ [0, 1] is an output drawn from the RNG. Advance the clock time to t t + x, and increment A(t) by 1 for t ≥ t .

2. If I(t) = n ≥ 1, the next event can be any of the three types. The time-interval until the next event is determined by x = log u/γ(n)T, where u ∈ [0, 1] as defined above. Take another RNG output u ∈ [0, 1]. Classify the event as type-(A), if 0 ≤ u ≤ P A : classify it as a type-(B) event, if P A < u ≤ P A + P B ; and classify it as type-(R), if P A + P B < u ≤ 1 Advance the simulation clock to t t + x, and increment the counter A(t), B(t) or R(t) by one depending on the classified result is type (A), (B) or (R).

3. Go back to step 1 or 2, depending I(t ) = 0 or I(t ) ≥ 1, and use t as a new clock time, and repeat the above steps.

As we discussed in Part I, the ratio of the standard deviation of a given probability distribution to its mean, called the coefficient of variation (CV), sometimes serve as a simple and good measure 4 The random number generator used today in MATLAB, Python, Pokemon, and others is what is known as Mersenne Twisted GFSR (general feedback shift register) sequence, developed by Makoto Matsumoto and Takuji Nishimoto (see "Mersenne Twister: a 623-dimensionally equidistributed uniform pseudo-random number generator," ACM Transactions on Modelling and Computer Simulation, January 1998). This algorithm generates a pseudorandom sequence of the period 2 19937 − 1, and the implementation is abbreviated as MT19937. It uses a Mersenne prime number M , of the form M = 2 n − 1, for some integer n. Sometimes n is restricted to a prime number, as well. 5 The property (a) has to do with the memoryless property of the variable X which has an exponential distribution. Assume that Y time units have elapsed since the last arrival. The time until the next arrival R X − Y has the same exponential distribution as X, regardless of Y . The property (b) means that when m independent Poisson processes with rate λ k , k = 1, 2, · · · , m are merged, the resulting stream is another Poisson process with rate m k=1 λ k . The property (c) is an opposite of the property (b). If a Poisson process with rate λ is split into m sub-streams, by assigning each arrival independently into the kth sub-stream with probability p k , where m k=1 p k = 1. Then the sub-streams are independent Poisson processes with rate p k λ's. For proofs of these properties, see e.g. [8] , pp. 403-405). 

97.5% 84% E[I(t)] 50% 16% 2.5% Figure 1 : Percentile curves of P n (t) of I(t)for 

97.5% 84% E[I(t)] 50% 16% 2.5% Figure 2 : Semi-log plot of the percentile curves

of the dispersion of the distribution. We obtained the CV of the distribution of I(t) in the BDI model, as given in (74) in Part I. If the distribution of a random variable X is close to a normal (or Gaussian) distribution, it is well known that realizations of X fall within the mean ±σ X with the probability 68.2%. Similarly, the mean ±2σ X provides a 95.4% confidence level, and the mean±3σ X gives 98% confidence level. But this rule of thumb does not apply to the NB(r, q) with small r, because the distribution is highly skewed, far from being symmetrical around the mean. A more accurate and reliable way is to compare simulation plots with the percentile curves of the NB(r, q). To that end, we calculate the cumulative distribution function (CDF) by

where P n (t) is the PMF of (16) and x is the largest integer, not exceeding x. Figure 1 shows the curves where the CDF take 0.975, and 0.025 (in red dash curves), 0.84 and 0.16 (in green dash) and 0.5 (in blue dash). The stochastic mean E[I(t)] I(t) is shown in black solid curve. We can expect that if we conduct many simulation experiments, about 68% of the sample paths will fall within the region between the green dashed curves, and about 95% of the sample paths should fall within the range given by the two red dashed curves. Roughly one half of the simulation curves should be above the blue dashed curve and the other half should be below this curve. Note that the stochastic mean curve is appreciably above the median (50%) curve: nearly by a factor of two.

Note that the top curve (the upper half of the 95% confidence interval) climbs up to as large as 100,000 by the 50th day, while it is only 1,800 on the 30th day. This is the power of the exponential growth, with which we are all familiar now by observing how rapidly the COVID-19's infections have grown in many parts of the world. This may not necessarily be a typical situation, because one half of the simulation runs, on average, should be above the "median curve", i.e., 50 percentile curve (shown in blue dash). Five of the six runs are within the 68% confidence interval, (between the two green dashed curves).

The most important observation to make here is the enormous disparity among the six simulations. Run 4 (in cyan) has as many as 1.22 × 10 5 infected people at t = 500 [days], whereas Run 2 (in red) has only ≈ 3 × 10 3 . Their ratio, therefore, is more than a factor of 40. If we run the simulator many times, for instance, the ratio of the numbers of the infected of the worst (e.g.,97% level) vs. the luckiest (e.g., 2.5% level) cases can be as large as 1,000 (see [9] , slide #33.).

If we plot the same simulation runs using the semi-log scale, as shown in Figures 6, we can see more clearly their behaviors in the initial rising phase (i.e., for small t) of these simulation runs of the stochastic process I(t). From the plots in the semi-log scale, we can see that once I(t) has reached the level of ≈ 100, the infected number in each run grows in a more or less deterministic fashion with a slope of 0.087 6 . This is because once the I(t) has reached ≈ 100, then a "(weak) law of large numbers" (As for detailed discussion on weak vs. strong law of large numbers, see e.g. [8] , pp.298-300.) will set in, because I(t) is the sum of many statistically independent and identical infection processes, and exhibits a more stable and predictable behavior than in the initial phase. In other words, the large variance among different sample paths of the BDI process is due to more erratic and unpredictable behavior in the early phase of the stochastic process I(t). In the initial phase the random arrivals of the infected from outside, as well as the fluctuations in the internal infection (a branching process) and the recovery/death, all contribute to the stochastic behavior of I(t) because I(t) is still small. As I(t) grows the effect of external arrivals will become negligible as long as r = ν λ is small, i.e., the order of unity or less. As remarked in the present author's keynote speech [9] , it may be worth noting that the BDI process model can be also used to explain the enormous disparity we observe in the wealth among different individuals of similar income levels, similar expenditure, and similar intelligence and knowledge in investment, thus they have a more or less similar investment portfolio in their initial phase. After 30 or 40 years after beginning their careers after schooling, some may become multi-millionaires (or even billionaires), whereas some end up with a life of modest means, if not living hand to mouth. The present author believes that this type of disparity in wealth distribution can be explained by applying the BDI process model. 7 . Figure 8 show the arrival patterns. The Poisson distribution of mean λ has the variance equal to the mean: σ 2 = λ; and for large λ, its CDF (cumulative distribution function) converges to that of the normal (or Gaussian distribution); Thus, 68 percent and 95 percent confidence levels could be well approximated by the λ ± √ λ, and the λ ± 2 √ λ, respectively. 6 log 10 e at = log 10 e × 0.2t = 0.087t. 7 In this analogy between the infection process and the individual's wealth growth, we interpret the infection rate λ as the ROI (return on investment), and the recovery/death rate µ as the rate of expenditures. The rate of arrivals of the infected from outside ν can be translated into the amount of additional money that can be put into investment. The parameter a = λ − µ is equivalent to the net growth rate of the wealth. Given that all these parameters are almost identical among different individuals, the disparity in their wealth growth depends, to a large extent, on their luckiness or unluckiness in the early phase of their investments, because the growth or decrease of wealth will be dictated by randomness in success or failure of investment, unexpected large expenditure or loss, and the availability of new fund for investment. Once the wealth exceeds a certain level, for instance, a million dollars, the future growth is more or less predictable from their investment strategy and expenditure because failures in some investments will average out with successes in other investments. Unavailability of new fund for investment will not affect much compared with the early phase when the wealth is small. The net asset growth will behave more or less predictably as a deterministic model can tell. To the best of the author's knowledge the applicability of the BDI process model to explain the disparity in wealth distribution seems a novel approach. The coefficient of variation of the Poisson process is given by

which converges to zero as t → ∞. Thus, the arrival process A(t) behaves much more predictably than the process I(t), whose coefficient of variation remains on the order of unity at all t.

c I(t) = 1

where r = ν λ and β(t) ≤ 1 is defined in (12) . As t → ∞, β(t) → 1, thus, the coefficient of deviation (CV) converges to λ ν , which, in our running example is 0.3/0.2 ≈ 1.225 as we discussed in Part I, p. 19, Example 3.

Although we were able to obtain the time-dependent probability mass function (PMF) of I(t), it seems rather difficult to obtain closed form expressions of PMFs for the processes B(t) and R(t), although we can find their PGFs. It appears that we need to be content with approximate PMFs of B(t) and R(t). The technique we will explore is the saddle-point integration method, which Bernhard Riemann (1826-1866) pioneered in his pursuit of the famous 1859 conjecture, known as the Riemann Hypothesis 8 . We will report our full discussion of the processes B(t) and R(t) in Part V [10] .

In the present section we only show the results of the above six simulation runs, together with the stochastic means B(t) and R(t), which we obtained in Part I, viz:

As we can see, the process B(t) and R(t) also exhibit enormous variations across the runs. They are positively correlated with each other and with the I(t) process, as well. Run 4 (magenta), Run 1 (blue) and Run 3 (yellow) are well above the mean curve, and Run 5 (green) and Run 2 (red) are below the mean in all the three processes B(t), R(t) and I(t). Run 6 (cyan) is just below the mean. The fact that these three processes behave in a similar fashion is quite expected, because both dB(t) dt and dR(t) dt are, on average, proportional 9 to I(t) at a given time. Referring to (2), A(t) is much smaller than I(t), B(t) and R(t), except for the initial period, thus we have an approximate formula

These large variations we observe in the processes B(t), R(t) and I(t) are all consequences of the built-in positive feedback loop inherent to the internal infection process B(t), which is a branching process 10 and gives rise to exponential growth exp(at).

In this section we discuss two additional topics for the benefits of the readers. One is how to derive the number of new infections for each day; the second is how to estimate the number of deaths. 9 New internal infections (excluding the new arrivals from the outside), occur at the rate of λI(t), and recoveries/removals/deaths occur at the rate of µI(t). Their expected values are related by the differential equations

as shown in Part I, page 11 (27) and (31). 10 See https://en.wikipedia.org/wiki/Branching_process and references therein.

The We take "day" as the time unit, as we do in our running examples throughout this report. We choose to interpret the interval t ∈ (0, 1) as the 1st day, then the number of newly infected persons reported on the tth day, denoted as N (t), should be computed as follows:

By substituting the B(t) and A(t) obtained in Part I, we readily find the stochastic mean or the expected value of I new [t]: ] by factor of 4 or 5; Run-6 (light blue) is close to the expectation. We observe enormous variations among the 6 runs. Since we have not yet obtained the PMFs for B(t), we cannot compute the percentile curves at this point. But from these simulation runs, the variations among different runs seem even more pronounced than we observed in I(t).

Note that for the short range 1 ≤ t ≤ 25, the vertical axis range is [0, 100], whereas for 1 ≤ t ≤ 50, the range is expanded by a factor of 150. This is consistent with the exponential growth of exp(at). From t 1 = 25 to t 2 = 25, it will grow by the factor of exp(a(t 2 − t 1 )) = exp(0.2 × 25) = e 5 = 148.4

The most important aspect of any model of an infectious disease should be how to estimate or predict the number of deaths. In the current COVID-19 epidemic, the case fatality rate is reportedly less than one percent for young people but will be much higher for 60 years or older, and those with comorbidity. In the present model, we assume a homogeneous population model, but the model can be generalized to multiple "classes" (or types), by assigning the model parameters ν c , λ c and µ c for different classes of infected population. As an illustration, let us assume fr= 2% as an overall fatality rate. Simulation of the death process D(t) can be done by randomly splitting the process R(t) with the specified fatality rate and produce D(t) as a sub-process of R(t). Figures 17 shows the result of the six runs. Figure 18 shows the first 25 days of the same D(t).

3 The BDI Process I(t) with the Initial Condition I 0 ≥ 1 By looking at the semi-log plots of simulation runs of the process I(t), some readers may wonder whether the huge variations among different runs may have to do with the fact that some simulation runs remain zero for a considerable period. This question can be answered by examining the process I(t) with the initial condition I 0 = 1. In referring to (9), the second term represents the PGF of the BD process without immigration (ν=0), which is often referred to as the simple birth-and-death process or as Feller-Arley (FA) process. 11

We rewrite (9) as G BDI:I 0 (z, t) = G BDI:0 (z, t)G F A:I 0 (z, t),

where G BDI:0 (z, t) is the PGF of the BDI process with the initial condition I 0 = 0, which we have studied thus far:

and

The product form expression (33) means that the BDI process with nonzero I 0 is the sum of the two independent processes, whose PGFs are given by (34) and (35):

By defining α(t) (see [13] [14] ) by

where β(t) was earlier defined by (12), we can rewrite (35) as

By taking the natural logarithm, differentiating w.r.t. z, and setting z = 1, we find the expected value of this random process:

Similarly, we find

which leads to the following expression 12 for the coefficient of variation (CV) of the FA process, when a > 0:

For the case of our running example with λ = 0.3 and µ = 0.1, the RHS of the above becomes 2/I 0 . When I 0 = 1, and the CV c F A:1 (t) ≈ 1.414 for all t, which is somewhat larger than the CV of the BDI process with I 0 = 0, which is c BDI:0 (t) ≈ r −1/2 = 1.225 (See Part I, p. 19 Example 2). 

We further analyze the above result depending (i) a > 0, (ii) a < 0, or (iii) a = 0.

(i) When a = λ − µ > 0: We find

where

is called the basic reproduction number in epidemiology, as defined in (32) of Part I. The distribution form of (43) is very similar to (16): it is a geometric distribution of the form ∝ β n (t), with β(t) ≈ 1, and has a long tail similar to Figures 7-12 in Part I. 12 The expression given in our earlier version [1] was incorrect.

(ii) When a < 0: We have lim t→∞ α(t) = 1, and lim

(iii) When a = 0, (i.e., µ = λ): We rearrange (37) as

Divide both the denominator and the nominator by a, and let a → 0. Using the familiar formula lim a→0 (e at − 1)/a = t, we find

Rearrange the PGF (35) in a similar fashion, and we find

from which we obtain

which could have been directly obtained from (39) and (40). The PMFs of (42) and (43) will be further simplified, when µ = λ:

Referring to (16), the PMF of the BDI process with I 0 = 1 is given by the convolution of the two PMFs, for which we have found closed form expressions:

In Figure 19 is the PMF P n (t) of the I(t) given in (16) at t = 10 for our running example, i.e. a = 0.2 and r = 0.667 (see Part I, P. 18, Figure 9 ). Figure 20 is the PMF (42) and (43). Figure 22 is the PDF of I(t) with the initial condition I 0 = 1 given by (16). We computed this by computing the convolution of the above two PDFs, i.e., (54). Figure 22 shows the CDFs of I(t) at t = 0, with I 0 = 0 (in blue) and with I 0 = 1 (in magenta). For the CDF F X (x) of a non-negative random variable X, you can show the following simple formula: 13

Thus, the white area should be equal to I(t) = ν a (e at − 1) = e 2 − 1 = 7.39 − 1 = 6.39 at t = 10. Thus, the blue bar area behind the magenta bars should be equal to I 0 e at = e 2 = 7.39, an increase in E[I(t)] at t = 10 due to the presence of an infected person at t = 0. An alternative and computationally more efficient way to compute the P BDI:1 n will be presented below. We use alternative representation of the PGF (35) as follow

where

Then we can have an alternative product representation to (33), as follows.

We may term the second term in the above product representation as the PGF of a generalized binomial distribution in the sense that the coefficients of z n terms are not necessarily non-negative. For I 0 = 1, this term is extremely simple, having only two terms: P 0 (t) = 1 − p(t) and P 1 (t) = p(t), which makes the convolution extremely simple: Fig 23 and Figure 24 show the PDFs of the above two PGFs. Their convolution results in the PMF as shown in Figure 21 . 

The percentile curves for the PMF P n (t) of I(t) with the initial condition I(0) = 1 are plotted together with the mean I(t) E[I(t)] are shown in Figure 25 and its semi-log plots are given in Figure 26 .

In Figure 27 we show six simulation runs with the initial condition I(0) = 1 and their semi-log plots in Figure 28 . The large variances among the sample paths still persist, thus we do not see any fundamental differences from the case where the simulations start with the initial condition I(0) = 0. This conclusion is not unexpected, if we go back to Equation (36):

I BDI:I 0 (t) = I BDI:0 (t) + I F A:I 0 (t).,

in which we set I 0 = 1. Since the two processes are statistically independent, the sum of their variances is the variance of the summed process:

where the first term on RHS is from Part I (53)

and the second term is from (40)

which leads to The expected value of the process I(t) with I(0) = 1 is from (39)

Thus, the CV of the the process I(t) remains constant for all t whether not the initial condition is zero or nonzero

which is even larger than the CV of I(t) with I(0) = 0, which is 1.225.

In the present report, we have presented simulation results by implementing a simulator of our stochastic model of an epidemic disease analyzed in Part I [2] . The simulated infection process I(t) indeed exhibits huge variations among different simulation runs. In this report, however, we presented only six consecutive runs in the interest of space. By presenting results of many more runs we could certainly show an even larger disparity between the worst and the best scenarios of a given BDI process I(t). By showing only several runs, however, we can sufficiently demonstrate the great value of a stochastic model in that a deterministic model is more often than not far from any sample path. In the same token, a single sample path of I(t), which we can observe in a real situation, can contain limited information about the ensemble of the process I(t).

In the current pandemic of COVID-19, most experts and policy makers in Japan seem looking for all possible factors that might help them understand why they see such enormous disparities between Japan and many countries in the world in terms of the magnitude of infected populations and death tolls. Our analysis and simulation shows that the epidemic process which is a branching process and is driven by inherently positive feedback can have enormous variations in the initial build-up phase, by mere luck or probabilistic chances, sometimes by factor of 100 or more, among environments with identical conditions, i.e., environments with the same λ, µ and ν parameters in our BDI process based model. In both analysis and simulations we have so far assumed a time-homogeneous model, whether the model parameters do not change in time. In Part III-A[3], we will report on our analysis of time-nonhomogeneous models, whereby we allow the model parameters λ(t) and µ(t) to be arbitrary functions of time. This generalization will help us better understand, for instance, how a change in social distancing, availability of effective vaccines and/or an improvement/degradation in medical treatment will affect the infection process. In Part III-B [4] we will report simulation results to help the reader better understand the significance of the analysis of the time-nonhomogeneous models and augment the analysis, which is limited in finding a closed form solution for a full-fledged BDI process.

In Part IV [5] , we plan to develop a statistical theory to estimate the model parameters from real data of COVID-19 epidemics and demonstrate how our stochastic model can be used in predicting the future behavior of an infectious disease, given its past and present value. Part V [10] will be devoted to use of saddle-point integration technique to approximately characterize the internal infection process B(t) and the recovery process R(t), which seem to defy an exact solution unlike the process I(t) for which we have an exact time-dependent probability distribution function.

Our BDI process based model has a fundamental advantage over nonlinear models such as SIR model and other models dominantly used in the epidemiology community in that our model should be readily extensible to more complex and realistic situations where different variants of the infectious diseases may coexist, and different classes of susceptible populations (e.g., aged population, those with comorbidity, etc.). These situations can be modeled by introducing a set of model parameters λ r , µ r and ν r , where r = (v, c) represent the variant v and the class c. ,

Stochastic Modeling of an Infectious Disease: Part II: Simulation Experiments and Verification of the Analysis

Stochastic Modeling of an Infectious Disease: Part I: Understand the Negative Binomial Distriution and Predict an Epidemic More Reliably

Stochastic Modeling of an Infectious Disease: Part III-A: Analysis of Time-Nonhomogeneous Models

Stochastic Modeling of an Infectious Disease: Part III-B: Simulation of Time-Nonhomogenous Models

Stochastic Modeling of an Infectious Disease: Part IV: Estimation of Model Parameters from Real Data, and Validation of Our Model

Modeling and Analysis: An Introduction to System Performance Evaluation Methodology

System Modeling and Analysis: Foundations for System Performance Evaluation

Probability, Random Processes, and Statistical Analysis

Stochastic Modeling of an Infectious Disease: Keynote speech at ITC-32

Stochastic Modeling of an Infectious Disease: Part V: Approximate Analysis of the Internal Infection Process B(t) and the Recovery Process R(t), based on Saddle-point Integration

Die Grundlagen der Volterraschen Theorie des Kampfes ums Dasein in wahrscheinlichkeitstheoritischer Behandlung

On the theory of stochastic processes and their application to the theory of cosmic radiation

The Elements of Stochastic Processes With Applications to the Natural Sciences

Lecture Note: Birth-and-Death Process and Its Application

Acknowledgments I thank Prof. Brian L. Mark of George Mason University for his valuable suggestions and help in the MATLAB simulation. Were it not for his kind help, I could not have completed this rather laborious study.