The spread of the corona virus and its infection as covid-19 has become a central issue in all of our lives. Curves showing the accumulated infections as well as the infection rate occur in the news every day. Most of us are 'social distancing' in order to 'flatten the curve' and we are receiving a variety of advice about how long it will take before it is safe to return to our normal activities. Consequently, I thought it would be both interesting and useful to calculate these curves using a simple model for the infection. It would be interesting just to see if I could actually do the calculation, and it would be useful if the solution made sense and gave me an appreciation of how the spreading takes place. The idea is not to perform a fully realistic calculation, but to learn what the relevant parameters are and to gain some insight about how they affect the spread of the disease. This paper is a summary of what I have learned.
Let us regard the affected population as a gas of colliding 'molecules' of a given number, N 0 , not nearly as big as Avogadro's number, N 0 = 6 × 10 23 atoms/mole, but ∼ 10 3 − 10 8 , comparable to populations of various regions of the country and the world. So we imagine that a population is like a gas of N 0 molecules, colliding with each other in a region of constant volume. Presumably, this gas will have characteristics of a real gas like temperature, pressure, and mean free path that can be determined. However, for the moment, we will bypass this interesting question and move forward to calculate the spread of the infection.
We assume the problem starts with N 0 − 1 molecules that are uninfected and one molecule that is both infected and contagious, in the sense that it will infect another molecule when it collides with that molecule. For simplicity, we suppose that there is an average time, τ s , between collisions, which we will use as the step time in our calculation. Referring to the left panel of Figure 1 , we begin at step 0 with the index j = 0 when a contagious (red) molecule is inserted into the gas of N 0 − 1 uninfected (blue) molecules. By the next step in the right panel of Figure 1 , the red molecule has collided with a blue one, making it red, so that two red molecules can go on to infect two other blue molecules in the next collision.
An important aspect of this problem is the finite time that a red molecule remains contagious. For the corona virus, we have heard that the lifetime is on the order of 2 weeks, the time that some potentially infected people have been told to quarantine. However, we will delay our treatment of this effect until later in this paper when we will have gained some insight and experience by solving the simpler problem for a permanently contagious virus. j = 0 j = 1 A gas of uninfected (blue) molecules with a contagious (red) molecule inserted at the initial time j = 0. (right) The same gas after one collision time when the red molecule has encountered a blue molecule and infected it (i.e. changed it to red).
Let us assume that the problem has advanced to the j−1 step at which point the number of infected molecules is represented by the symbol, n (j−1) inf , where the exponent, j − 1, is a tag and not a power to which n inf is to be raised. The question is then what is the number of infected molecules in the next step, as represented by the symbol n (j) inf ? For the first few steps with j ∼ 1, the probability, p j−1 , of colliding with an uninfected molecule will be close to 1 so that n (j) inf would be 2n (j−1) inf (i.e. just twice the present number of infected molecules before the collision). However, for larger values of j, we must consider the case of p j−1 = 1.
The probability that a red molecule will collide with another infected molecule is just (n (j−1) inf − 1)/(N 0 − 1), so the probability that it will collide with an uninfected blue molecule is
Consequently, n
and therefore n (j)
Now, for simplicity, we normalize the number of molecules to the total number, N 0 , in the gas, and define ν j = n (j) inf /N 0 as the normalized quantity. Thus, ν 0 = 1/N 0 , and
Next, in Figure 2 , we plot the growth curve, ν j , and the growth rate, ∆ν j = ν j − ν j−1 , versus the number of steps (j) or collision times (jτ s ) that have elapsed since the start of the infection. The difference equation Eq. which is one collision before the growth rate in the right panel reaches its peak value of 0.25 at j max ≈ 17. At this point, 1/4 of the total population is being infected per collision time. The full width at half maximum is approximately 4 collision times, which means that the area under the curve is about 1, consistent with the final height of the infection curve in the left panel. Also, note that each profile is asymmetric. In the left panel, the trip from the point of maximum slope to the upper level is shorter than the trip from the lower level to the point of maximum slope. Likewise, in the right panel, the trip from the peak to the final base level is shorter than the trip from the initial base level to the peak.
Although the shape of the growth-rate profile does not vary with N 0 , its location, j max , does vary with N 0 , as shown by the measurements that are plotted as red dots in Figure 3 . The least-squares best fit to these measurements is shown by the dashed red line, whose This red line is nearly coincident with the dashed blue line, which represents the relation j max = log 2 N 0 = log 10 N 0 log 10 2 ≈ 3.322 log 10 N 0 , (3.6) -8which is the number of two-factor collisions that it takes to reach the total population, N 0 .
It is instructive to see if we can find an analytical solution to Eq. (3.4). Let's begin with the growth equation in the form
Initially, ν j−1 << 1, so that ν j ≈ 2ν j−1 = 2 2 ν j−2 = · · · = 2 n ν j−n . Setting j = n, we obtain the solution ν n = 2 n ν 0 , corresponding to an exponentially increasing number of infections. So it is the second term that limits the exponential growth and stops the infection.
We can use a similar approach when the second term is included. In this case, we change variables by introducing µ j = 1 − ν j and obtain
Consequently, µ j = µ 2 j−1 = µ 4 j−2 = µ 8 j−3 = ··· = µ 2 n j−n . Setting j = n, we obtain µ n = µ 2 n 0 . Finally, changing back to ν j , we have
as the solution to Eq.(3.7).
Next, noting that 1 − ν 0 ≈ e −ν 0 , we can change Eq.(3.9) to ν j ≈ 1 − e −ν 0 2 j = 1 − 2 −(ν 0 /ln2)2 j = 1 − 2 −2 (j−k) , (3.10)
where (ν 0 /ln2) = 2 −k . Because ν 0 = 1/N 0 , it follows that k = log 2 (N 0 ln2), which is 1 unit less than j max , the location of the peak of the growth rate, ∆ν j = ν j − ν j−1 , as we shall see next.
We begin by using Eq.(3.7) to evaluate ∆ν j = ν j − ν j−1 . The result is
The quantity, ν j−1 (1 − ν j−1 ), has its maximum value of 1/4 when ν j−1 = 1/2. Thus, ∆ν j has its maximum value when ν j−1 = 2 −2 j−1−k = 1/2. This happens when j − 1 − k = 0, and therefore j max = j = k + 1 = 1 + log 2 (N 0 ln2), (3.12) which is j max = log 2 N 0 + 0.471 when the numerical value of 1 + log 2 (ln2) is used. Also, because ν j−1 = (1 − ν 0 ) 2 j−1 ≈ 2 −2 (j−jmax ) , it follows that ∆ν j = 2 −2 (j−jmax ) 1 − 2 −2 (j−jmax ) . Summarizing the results for both ν j and ∆ν j , we have
As we saw in Figure 2 , these solutions for ν j and ∆ν j agree very well with the data points that were obtained from the recursive equations, Eq.(3.4) (or (3.7)) and Eq.(3.11). The point at j = 17 occurred very slightly to the left of the peak at j max = 17.1. Also, we previously noted in the right panel of Figure 2 that the full-width at half maximum (w) of the growth-rate peak is approximately 4 collision times. Now, with Eq.(3.10), we can calculate this width precisely by setting ν j−1 (1 − ν j−1 ) = (1/2)∆ν max = 1/8. This is a quadratic equation whose roots are ν j−1 = (1/2)(1 ± 1/ √ 2). Setting them equal to 2 −2 (j−k) , we obtain 2 j−k = 1 − log 2 (1 ± 1/ √ 2), whose numerical values are 0.2284 and 2.7716. The corresponding values of j − k are −2.130 and 1.471, and their difference (the full-width at half maximum) is w = 3.601, about 10% less than our original estimate from Figure 2 .
It is interesting to compare the width and height of this growth-rate curve with the total area under it (which we expect to be 1 because j=25 j=0 ∆ν j = ν 25 − ν 0 = 1 − 1/N 0 ≈ 1). For our calculated full-width at half maximum, the product of height and width is 0.25×3.601 = 0.900, which is slightly less than 1. On the other hand, if we had used the full-width at 1/e maximum instead of the full-width at half maximum, we would have obtained 0.25×4.398 = 1.100, which is slightly greater than 1. It is an interesting coincidence that the average of these two widths is almost exactly 4.000, which gives an area of 1.000. Now, let us return to the more realistic case in which the red molecules are contagious for a finite time.
In this section, we shall assume that the infected red molecules do not remain contagious indefinitely, but lose their ability to infect other molecules after a finite time, c. We suppose that c > 1, so that the red molecules remain contagious for at least one collision time. Otherwise, the molecules would lose their ability to infect before they encounter uninfected blue molecules in the next collision, and the epidemic would be over. So, we suppose that after a red molecule suffers c collisions, it loses is ability to infect others. At this point, we change its color from red to green.
In the real world, a molecule could become green either because it recovered from the spreading disease or because it did not recover and died. In the latter case, the green molecule ought to be removed from the gas. For the current corona virus, this applies to about 3-5% of the green molecules (or less depending on how many of the infections were not reported). Consequently, in this paper we will neglect this effect and keep all of the green molecules in the gas. As we shall see later, these green molecules provide a kind of 'shielding effect' by reducing the probability that the red molecules will come in contact with the remaining blue molecules during the next collision. Figure 4 shows the progress for c = 2. In the left panel, one collision has taken place, and the original red molecule from j = 0 has infected another molecule at j = 1. These two red molecules continue on to infect two more molecules by j = 2 (right panel), but by this j = 1 j = 2 time, the original red molecule has lost its contagiousness and has become green. Like the three red molecules, the green molecule will go on to collide with other blue molecules, but it will not infect them or change their color.
Proceeding as we did in the previous section, the probability that a red molecule will collide with a blue molecule is still given by Eq. (3.1). However, because some of the infected molecules are no longer red (i.e. contagious), Eq. (3.2) must be replaced by
where n (j−1) c is the number of contagious molecules at the time j − 1. Now, we need to find an equation for the number of contagious molecules.
To get a feel for this process, we refer to Figure 5 , which is a sketch of the evolving number of red and green molecules as a function of the collision time, j, assuming that there are no collisions between red molecules during this time. At the left end of the array, we begin at j = 0 with R 0 , which is a symbol for a single red molecule of 'age' 0. After the first collision at j = 1, there are two red molecules -the original red molecule, R 1 , whose age is now 1, and the newly infected molecule, R 0 . One step later, there are 4 infected molecules -R 0 and R 1 , derived from R 0 at time j = 1, plus R 0 and G, derived from R 1 at time j = 1. Instead of calling the the originally infected molecule by the name R 2 , we realize that this molecule is not contagious anymore, and we represent it by the letter G. As time passes to j = 3 and j = 4, the pattern expands, with each R 0 always producing an R 0 and an R 1 , with each R 1 producing an R 0 and a G, and with each G always remaining a G.
From the evolving pattern in Figure 5 with c = 2, we can see that the number of contagious molecules at time j − 1 is equal to the total number of infected molecules (contagious plus noncontagious) at that time minus the total number of infected molecules 2 steps earlier.
More generally, we can appreciate that all of the contagious molecules (R's) at any time, j − 1, will lose their contagiousness (and become G's) after an additional time, c, so that
(4.2) even when there are collisions between red molecules. Combining Eqs. (3.1), (4.1), and (4.2), and normalizing to the total population, N 0 , we have
which reduces to Eq. (3.7) when ν j−1−c = 0. Likewise, the difference form of Eq. (4.3) is
which reduces to Eq. (3.11) when ν j−1−c = 0.
Our next step is to plot ν j and ∆ν j versus j for various values of the parameter, c. To do this, we set ν 0 = 1/N 0 and use Eqs. (4.3) and (4.4) for j ≥ 1. Also, we add the condition that ν j = 0 for j < 0. This means that ν j−1−c will not contribute until j − 1 = c.
Finally, because the influence of the finite contagious time becomes greatest as that time approaches the collision time (i.e. as c → 1), we need to find a way to make these equations apply for non-integral values of c between 1 and 2. For this purpose, we used a linear interpolation between values of c in the range [1, 2], or more generally in the range [c 0 , c 0 + 1], where c 0 is an integer. Consequently, we replaced ν j−1−c by the expression λν j−1−c 0 −1 + (1 − λ)ν j−1−c 0 , where λ lies in the range [0, 1] and c = c 0 + λ.
Figures 6 and 7 show ν j and ∆ν j plotted versus j for parametric values of c from 1.25 to 10.0. I chose 10 as the largest contagious time because the curves for c = 10 and a much larger value of c like c = 100 were indistinguishable. Thus, once the contagious time reaches 10 collision times, the result is the same as if the red molecules remain contagious for the duration of the problem. (The curves were almost indistinguishable for c = 6, except that the peak height for the differential curve was 0.24 instead of 0.25.) Consequently, the curves with c = 10 are an accurate representation of the result when the red molecules remain contagious for the duration of the problem (i.e. the value of ν j eventually reaches 1, and the height of the peak of ∆ν j is 0.25 at the location where j max ≈ 17, which is approximately log 2 N 0 = 16.6.)
In retrospect, it is obvious why the growth curve with c = 10 is indistinguishable from the corresponding curve for permanently contagious molecules (with c = ∞). For c = 10, the lag is so long that the plot of ν j−1 finished its rise from 0 to 1 before the lagged quantity, ν j−1−c , became appreciable (ν 9 = 0.005). Equivalently, there are 10 points between ν j ≈ 0 and ν j ≈ 1. For a smaller value of c, the lag is shorter, and ν j−1−c becomes appreciable before ν j reaches the end of the rising portion of its curve. In addition, there are more points between the beginning and end of the rising portion of the curve. This produces less spacing between the points of ν j in Figure 6 and lower peaks in the plots of ∆ν j in Figure 7 .
For c = 10, the fraction of infected molecules reaches 1 at j ≈ 20 and remains there forever. This means that everyone in the population eventually becomes infected. However, as c decreases, the final fraction of infected molecules decreases, reaching about 0.38 when c = 1.25 collision times. If the contagious time were about 14 days, c = 1.25 would correspond to a collision time of 11.2 days. This means that to reduce the number of infected molecules to 38% of the total population, one would have to increase the average time between collisions to about 11 days. Such 'social distancing' prevents 62% of the population from becoming infected and shows the importance of making the time between collisions as close as possible to the 14-day contagious time. But even for 11 days, the end of the epidemic would be delayed appreciably from j ∼ 20 to j ∼ 70, which is roughly from 1 year to 3 years if the collision time is 2 weeks. Figure 7 shows that the infection rate weakens, broadens, and shifts to later times as the contagious time, c, approaches the collision time. I think this is what the media reports for the corona virus are referring to when they use the term 'flattening the curve'. Corresponding plots for other values of N 0 show essentially the same profiles, shifted in time. In the next section, we shall analyze the dependences of these profiles on N 0 and c.
In the previous section, we saw that the growth rate gradually rises to a peak ∼ 0.25N 0 τ s −1 when the contagious time, c, is large compared to 1, and that the peak height becomes smaller as c approaches 1. In this section, we will determine the location, height, and width of this peak as a function of N 0 and c. Also, we will determine the dependence of the final value of ν j on these parameters.
We can obtain a clue to the relation between c and ν f by examining the behavior of ν j−1 in the region where the growth curve has its maximum slope. In this region, the slope is fairly constant so that it will have approximately the same value at several neighboring points. In particular, we would expect that
if c is close to 1. In this case, Eq. (4.4) reduces to
which would be an identity if c(1 − ν j−1 ) = 1. Equivalently,
This point of maximum slope lies somewhere between the starting and ending values ν 0 and ν f , perhaps becoming centered between them if the growth-rate curves become symmetric when c approaches 1 (as Figure 7 suggests). In this case, we would expect that
where α is a number on the order of 2.
Although Eq.(5.4) seems to reproduce the calculated values of ν f fairly well as c approaches 1, it does not satisfy the requirement that ν f → 1 as c → ∞. To satisfy this requirement and still retain the behavior for c ≈ 1, we can use an exponential of the form 1 − e −α(c−1) . This expression approaches 1 as c − 1 → ∞, and it approaches α(c − 1) as c − 1 → 0. In fact, if we retain only the first-order terms in c − 1, we obtain the equality: where the final value of ν f occurs, but it does not determine ν f itself. In other words, we would obtain the same plot of ν f versus c for other values of N 0 . The solid red curve is the least-squares-best-fit to these points using the relation ν f = 1 − e −1.890(c−1) , (5.6) and the dashed blue curve is the approximation given by Eq. (5.4) with α = 1.890.
As we can see, Eq. (5.4), represented by the dashed blue line, fits the data very well for c up to 1.5, but the exponential relation given by Eq. (5.6) is required to extend the fit to larger values of c. Consequently, Eq. (5.6) provides the desired relation between the amount of social distancing and the fraction of molecules that eventually become infected (and therefore the fraction that ultimately escapes the epidemic). Moreover, for social distancing with c < 1.5, a rough estimate of the final fraction of infected molecules is just ν f ∼ 2 (1 − 1/c).
In the previous section, we found that Eq. (5.6) provides an accurate expression for the final fraction of infected molecules, ν f , as a function of the parameter, c. Next, we look for analogous expressions for the height and width of the peak in the growth-rate curve.
We begin by referring to Figure 9 , which contains plots of the accumulated infection, ν j , and the infection rate, ∆ν j , as a function of the collision step number, j. The left panel shows the accumulated infection curve (in red) and a black line segment drawn tangent to it at the point of maximum slope. Actually, this point is really an extended region, which is why we were able to deduce that ν j−1 ≈ 1 − 1/c in the center of this region. (See the discussion leading to Eq. (5.3).) The ends of this line segment lie on the horizontal blue lines where ν j = ν f and ν j = ν 0 = 1/N 0 ≈ 0, and are separated by a horizontal distance, ∆j. Consequently, ν f /∆j is the maximum slope of the ν j curve, which in turn is equal to the peak value, ∆ν max , of the growth-rate curve in the right panel. (Here, we ignore the subtle difference between ∆ν j and dν j /dj because ∆ν
where h is the peak height, ∆ν max .
Next, in Figure 10 , we plot the numerical values of the ratio h/ν f versus the parameter, c, and fit the resulting points with two slightly different formulas given by Because ∆j is a measure of the width of the growth-rate curve, we refer to it as the width, w e , of that curve. Consequently, that width can be obtained by taking the reciprocal of the expressions in Eqs.(5.8a) and (5.8b) 1 . As we can see, the simpler formula given by Eq.(5.8a) ends up at 0.25, but misses some of the points at the bend of the curve. However, the more complicated formula given by Eq.(5.8b) passes through the points at the bend of the curve, but ends up slightly lower than 0.25. Both fits pass through all of the points where c ≤3.
The next step is to see how well these expressions combine with Eq.(5.6) for ν f to fit the calculated points for h. The result is shown in Figure 11 . Here, the product of the two fits reproduces the initial curvature in the plot of h versus c. Also, as expected, the better 1 Other measures of the width like the full-width at half maximum are slightly smaller than w e . For example, the width, w, would be given by w = ln(1 + √ 2)w e = 0.881w e if the fraction of infected molecules, ν j , were fit by a hyperbolic tangent of the form ν j = (ν f /2)[1 + tanh{γ(j − j max )}], and the width would be w = (2 2/3 − 1) 1/2 w e = 0.766w e if ν j were fit by a square-root expression of the form ν j = (ν f /2)[1 + γ(j − j max ){1 + γ 2 (j − j max ) 2 } −1/2 ], where γ = 2/w e . fit to 1/∆j in Figure 10 provides the better fit to h in Figure 11 . However, even the simpler expressions given by Eqs.(5.6) and (5.8a) provide an essentially perfect fit for c ≤ 3, which is the region of interest when a substantial amount of social distancing is maintained.
In summary, we can describe the growth curves of Figure 6 and growth-rate curves of Figure 7 in terms of the contagious lifetime parameter, c, using the approximate relations for ν f , 1/w e , and h:
(5.9c)
I have used plots like those shown in Figure 7 to measure the locations, j max , of the peaks for values of c ranging from 1.25 to 10.0, and values of log 10 N 0 ranging from 2 to 8. Figure 12 shows the measured values of j max plotted versus log 10 N 0 for given values of c. The dashed lines are the least-squares best fits to these points given by: 
−14.071 + 12.357 log 10 N 0 , c = 1.25, −5.154 + 7.439 log 10 N 0 , c = 1.5, −1.381 + 4.799 log 10 N 0 , c = 2, −0.156 + 3.796 log 10 N 0 , c = 3, +0.312 + 3.359 log 10 N 0 , c = 10.
(5.10) Also, the black line represents Eq.(3.6) (j max = log 2 N 0 ≈ 3.322 log 10 N 0 ) for permanently contagious molecules with c = ∞. Its slope (3.322) and intercept (0.000) are very close to the slope (3.359) and intercept (0.312) of the line with c = 10.
Based on this figure, there can be little doubt that j max = p c + q c log 10 N 0 for each value of the parameter, c, where p c and q c are the least squares best-fit parameters for the dashed lines that pass through these points. However, an even more interesting property of Figure 12 is that the spread of the j max values increases with log 10 N 0 , and that this spread tends to a very small minimum value near the point (1.5, 6) where the straight lines nearly intersect. This means that for a small population, the maximum growth rates occur at approximately the same time regardless of the value of c and the corresponding amount of social distancing. This contrasts with the result for a large population for which the location of the peaks depends much more strongly on the amount of social distancing. Thus, to achieve the same maximum growth rate, the larger population must wait a much longer time than a smaller population. As we shall see in subsection (5.3.2), this property can be used to deduce j max as a function of log 10 N 0 and c. But first in subsection (5.3.1), we will deduce an expression for j max using the linearized solution of the growth equation given in Appendix A.
Appendix A provides a way to estimate the values of p c and q c and thereby reproduce the j max expressions in Eq. (5.10). Referring to Figure 25 in Appendix A, we see that the linearized solutions give the initial behavior of ν j for each value of c. Because these linearized solutions lack the damping effect of the non-linear term, they eventually diverge from the nonlinear solutions whose curves bend over and approach final values of ν f . However, the plots of the linearized solutions in Figure 25 seem to cross the ν j = ν f threshold approximately where the non-linear solutions have their maximum slopes. Such a relation would provide a way of determining the location of those maximum slopes and therefore the locations of the peaks of the corresponding growth-rate curves of ∆ν j . For this purpose, we would simply equate the linearized solution, σ c ρ j c ν 0 , and the final value of the growth curve, ν f , and then solve for j. The first step gives
where ν f is obtained from Eq.(5.6), and the parameters σ c and ρ c can be obtained from Table  1 of Appendix A (or from Eq. (A9) for values of c in the range 1-2), or from Eqs. (C1a) and (C1b) of appendix C. The resulting value of j is j max − 1 because this analysis is based on the growth curves in Figure 6 and Figure 25 , whereas the peak value of the growth-rate curve comes from the difference, ∆ν j = ν j − ν j−1 shown in Figure 7 . However, because this method is based on a rough estimate, I would be surprised if the accuracy of the resulting value of j max were closer than 1 unit of collision time. This is comparable to the differences between ∆ν j and dν j /dj, which Figure 7 shows are often about 0.5 of a unit.
Solving Eq. (5.11) for j and setting j max = j + 1, we obtain
which is in the format p c +q c log 10 N 0 used in Eq. (5.10). Although, I have used logarithms to base 10 for comparison with the expressions in Eq. (5.10), because these logarithms occur as ratios, they can be expressed in any base, including base e, base 2, or even base ρ c , without affecting the value of j max .
Next, obtaining ν f from Eq. (5.6) and substituting the values of σ c and ρ c from Table  1 of Appendix A, we find that
934 + 12.239 log 10 N 0 , c = 1.25, −4.802 + 7.383 log 10 N 0 , c = 1.5, −0.668 + 4.785 log 10 N 0 , c = 2, +0.463 + 3.780 log 10 N 0 , c = 3, +1.000 + 3.322 log 10 N 0 , c = ∞.
(5.13) These values are comparable to the least-squares best-fit values in Eq. (5.10). The coefficients of log 10 N 0 are essentially the same. However, the constant terms in Eq. (5.13) are systematically a little larger than those in Eq. (5.10). For the larger values of c, this discrepancy is probably due to the intrinsic differences between the digital peak and the peak of the continuous curve. For the smaller values of c, and in particular for c = 1.25, the discrepancy may reflect the weaker and broader profiles whose peak locations are less well defined. Figure 13 shows the j max values calculated using the full non-linear equation (Eq. (4.4)) and the values obtained with the approximate method given by Eq. (5.12), both using the nominal value of log 10 N 0 = 5. The red dots are based on Eq. (4.4), whereas the blue dots and their interpolated curve are based on Eq. (5.12). The red and blue points all lie on the same dashed blue curve, indicating a clear dependence of j max on c. For c ≥ 1.4, the corresponding red and blue points of each pair are very close together, often occulting each other with virtually identical values. However, for c < 1.4, the curve steepens and the separations systematically increase as c approaches 1. At first, one might suppose that the overall agreement between the calculated values in Eq. (5.13) and the measured values in Eq.(5.10) supports our ansatz that the point of maximum slope occurs where the linearized solution for ν j equals the final value ν f . However, the agreement was better for q c than for p c , which contains the entire dependence on ν f . So the discrepancy in the values of p c probably reflects the eventual breakdown of this assumption, especially as c approaches 1 and the growth-rate curves become increasingly flat. Next, we will abandon this assumption, and see if the nearly common intersection of the lines in Figure 12 gives a point-slope formula for j max that is more accurate than the slope-intercept approach.
In this approach, we return to Figure 12 where the measured values of j max give linear fits that nearly intersect at a common point near (x 0 , y 0 ) = (1.5, 6). We begin by precisely determining the point of closest approach to the lines with c = 1.25, 1.5, 2, 3, and ∞, leaving out c = 10 because it is essentially the same as c = ∞.
We do this in the following way. We recognize that even though these five lines may not intersect in the same point, lines drawn parallel to them can be adjusted so that they do intersect at a common point. This can be done for many such common locations and our objective is to find the one that makes the root-mean-square distance between all those parallel lines a minimum. (Because one set of lines passes through a common point, (x 0 , y 0 ), this is equivalent to making the mean-squared distance from that point to the original lines a minimum.) Consequently, we derive an expression for the mean squared distance, D 2 (x 0 , y 0 ), between all of the parallel lines. Then we minimize D 2 (x 0 , y 0 ) with respect to x 0 and y 0 , and obtain two linear equations for x 0 and y 0 . When we solve those equations simultaneously, we obtain x 0 = 1.586 and y 0 = 6.153. Using these values to evaluate D 2 (x 0 , y 0 ) and taking the square root, we find that D rms = 0.054.
Next, we write the equation for the lines that pass through the point (x 0 , y 0 ) with the same slope, q c = 1/ log 10 ρ c , that we used in the previous subsection. This point-slope approach gives the result
where (x 0 , y 0 ) = (1.586, 6.153) and q c = 1/ log 10 ρ c . The quantity, ρ c , is given by Table 1 of Appendix A. It is the largest positive root of the equation r c+1 − 2r c + 1 = 0 when c is an integer greater than 1. However, when c is not an integer, it is necessary to interpolate the linearized difference equation to obtain an equation for r, as we did in Appendix A for c between 1 and 2. An approximate expression is given in Eq.(C1a) of appendix C. Figure 14 shows these lines plotted with the measured points from Figure 12 . The passage of these lines through the common point, (x 0 , y 0 ), is clearly visible. These lines fit the measured points almost as well as the best-fit lines in Figure 12 , verifying that j max is well represented by a line of slope 1/ log 10 ρ c , passing through a point, (x 0 , y 0 ), that does not depend on c. Eq.(5.15) gives numerical values for comparison with Eqs.(5.13) and (5.10). 
It is interesting to interpret the growth equation in terms of the fractions of red, green, and blue molecules. For j ≤ c, all of the infected molecules are still contagious (red) and the incremental increase of infections (normalized to the population, N 0 ) from one collision to the next is
Here, ν j−1 is the fraction of infected molecules at time j −1 (all of which are red) and 1−ν j−1 is the fraction of noninfected blue molecules, which is also the probability of colliding with a blue molecule. Consequently, the incremental increase in infections is just the product, ν j−1 (1 − ν j−1 ). However, for j ≥ c + 1, the molecules that were infected c steps earlier are no longer contagious. Consequently, at step j − 1, the fraction of infected, but noncontagious, green molecules is ν j−1−c , and the incremental increase becomes
Thus, after a time, c, the increase in the number of infections is still equal to the product of the numbers of red and blue molecules, but this increase is smaller than it would have been if the red molecules were permanently contagious and there were no green molecules. Eq.(6.2) indicates that the increase is smaller because the number of red molecules is reduced by the current number of green molecules.
It is also interesting to consider what happens toward the end of the epidemic when ∆ν j → 0. For permanently contagious molecules, Eq.(6.1) tells us that the factor (1 − ν j−1 ) → 0 and therefore the other factor ν j−1 → 1. In this case, the epidemic ends because all the blue molecules become infected, and now exist as contagious red molecules. However, if the infected molecules remain contagious for a limited time, c < 10, then all of the blue molecules do not become infected and the factor (1 − ν j−1 ) does not approach 0. Instead, the epidemic ends when other factor (ν j−1 − ν j−1−c ) → 0. This means that all of the infected molecules become green and there are no red molecules left to infect blue ones. As we have seen in Figures 6 and 8 and Eq.(5.6), by decreasing c (and thereby increasing the amount of social distancing), we cause a greater number of blue molecules to avoid infection by the time that the epidemic is over. Figure 15 shows the separate evolutions of red (ν j − ν j−c ) and green (ν j−c ) molecules as well as their sum, ν j (the combined total of all infections, whether still contagious or not). (We do not need to provide a separate plot for blue molecules because their number, 1 − ν j , is just 1 minus the the combined number of red and green molecules, ν j .) The dashed black curves show the familiar 3-part variation of ν j as a function of time in units of the step time, τ s . The change begins slowly for a while, then increases rapidly for a time on the order of w e = ∆j, and then levels out at a final fraction of infected molecules, whose value, ν f , depends on the value of c. Because the number of green molecules is given by the lagged quantity, ν j−c , it shows the same time dependence as the number of infected molecules, ν j , but delayed by c collision times.
The red molecules do not have this 3-part monotonic increase. Shown by the red curves in Figure 15 , the number of red molecules have a peaked profile, increasing to a maximum value and then falling back to zero at the end of the epidemic. This is what one would expect from their definition, ν j − ν j−c , which is the difference between two terms that eventually become equal.
We can learn more about the time variation of the red molecules by examining the plots in Figure 15 as a function of c. When c is much smaller than the width of the peak (which is essentially the rise time, w e = ∆j, of the black profile of total infections), we see that
Thus, for c << w e , the plots of the number of red molecules are within a factor of c of matching the plots of infection rate, ∆ν j , in Figure 7 . This near equality is shown in Figure 16 , which compares ν j − ν j−c and c∆ν j for c = 1.5 and 2 and log 10 N 0 = 5. In the left panel, the two curves for c = 1.5 are nearly identical. In the right panel where c = 2, the peak of the red curve is still nearly equal to the height of the black curve, but is now shifted to the right by about 0.5 of a collision time. However, the agreement breaks down as c becomes larger and the steep increase in the number of green molecules is appreciably delayed. For c ≥ 3, we return to Figure 15 and see that the rising segment of the red curve coincides with the rising segment of the black curve, until the infection rate is dominated by green molecules and the red curve reaches its peak. This leads to the concept of 'herd immunity' which we shall discuss next. 
An equivalent way of thinking about the relation between the numbers of red and green molecules is through the so-called 'herd immunity' concept that one often encounters in the literature (see, for example, Fine et al. (2011) ). Here, we think of the green molecules as shielding blue molecules from the remaining red ones, so that some of these red molecules cannot find blue molecules before they lose their contagiousness and turn green. The increased number of green molecules makes it even more difficult for the surviving red molecules to find blue ones, and the number of red molecules starts to decrease. We take the threshold for this turnaround to be the number of infected molecules (relative to the total population, N 0 ) when the number of red molecules has reached its peak. However, as Eq.(6.2) indicates, the remaining blue molecules are still subject to infection as long as there are red molecules in the population, and the epidemic is not over until all of those red molecules have turned green. Consequently, even after reaching herd immunity, an appreciable time remains before the red population is reduced to 0 and a blue molecule is safe from infection.
We can think about this process graphically with the help of Figure 15 . The number of red molecules is given by ν j − ν j−c , the difference between the black and green curves. This difference is maximum where ∆(ν j − ν j−c ) = 0, which means that ∆ν j = ∆ν j−c . Consequently, the red curve reaches its peak where the profiles of ∆ν j and ∆ν j−c intersect, which is somewhere between the peaks of those two profiles, depending on the value of c. Another way of interpreting the equality, ∆ν j = ∆ν j−c , is to say that the red curve reaches its peak at the value of j where the slopes of the black and green curves are equal. When c is small, these points lie slightly above and below the points of maximum slope on the black and green curves, corresponding to a relatively small peak in the red curve. As c increases, the peak height of the red curve increases, and the points of equal slope on the black and green curves move farther apart, so that ν j increases and lies closer to the 'knee' of the black curve and ν j−c decreases and lies closer to the 'toe' of the green curve. Eventually, for c = 10, ν j = 1 and ν j−c = 0, corresponding to a common slope of 0.
To describe this process more quantitatively, we refer to Figure 17 , which compares plots of the numbers of red molecules (ν j − ν j−c ) with corresponding plots of the total number of infections (ν j ) for values of c in the range 1.25 to 20. In the upper panel, we find the time that a particular plot of the number of red molecules reaches its peak, and then in the lower panel, we find the value of ν j for the corresponding plot of the total number of infections at that time. Thus, in the upper panel, the blue curve with c = 1.5 has a peak of strength 0.08 at j = 33 which corresponds to the value ν j = 0.38 in the blue profile in the lower panel. As c increases, the threshold value of ν j also increases, eventually approaching 1 as c becomes large. In the lower panel, the curves with c = 5 and c = 20 are almost indistinguishable, indicating that ν j is already very close to 1. This is consistent with our earlier observation that the growth is essentially the same for c = 10 as it is when c → ∞. Figure 18 shows a plot of these computed values of ν j , which we now refer to as the herd-immunity threshold, ν herd , versus the corresponding values of c. The solid red curve is the root-mean-square best fit to those points using the formula
Having learned what the fraction of total infections is when the corresponding fraction of contagious red molecules reaches its peak, we now ask when the fraction of red molecules actually reaches its peak. In particular, after reaching the time, j max , that the total number of infections (red plus green molecules) reaches its peak, how much longer must we wait until the number of red molecules reaches its peak? In Figure 16 , the plots of ν j − ν j−c and c∆ν j have already shown that the wait time is roughly 0.5 of a collision time when c = 2 and even less when c = 1.5. So the next step is to make additional comparisons and plot the results. Figure 19 shows the result of measuring the peak locations of ∆ν j and ν j − ν j−c . The measured lags are relatively small, never exceeding 2 collision times, which is the difference, -The time lag, j red − j max , between the peaks of the total number of infections, ν j , and the number of contagious molecules, ν j − ν j−c , plotted versus the contagious time, c. The solid red curve is the root-mean-square best fit to the calculated points (red dots) given by Eq.(6.5), and the dotted red line is its extension with a constant lag of 2 after the solid curve reaches its peak at c m = 8.81.
j red − j max , when c → ∞. (For this calculation, we return to Figure 2 and recognize that j red = 19 is the j-value when ν j reaches 1, as shown in the left panel, and that j max = 17 is the time that ∆ν j reaches its peak, as shown in the right panel.) Coming back to Figure 19 , we note that the data points are best fit by the quadratic function
where c m is the value of c where this quadratic function reaches its maximum value of 2 units of the step time, τ s . In this case, c m = 8.81. The numerical form of this equation shows that the lag starts as a linear function whose slope is nearly 1/2, but 'decelerates' very slightly as c becomes large. Taking the derivative of Eq.(6.5), we obtain the slope
which equals 0.51 when c = 1 and becomes 0 when c = c m . Thus, if c − 1 is small, then j red − j max ≈ 0.51(c − 1), which gives lags of 0.51 for c = 2 and 0.25 for c = 1.5, consistent with our measurements from Figure 16 .
It is interesting to compare the c-dependence of this lag with the corresponding formula obtained for profiles like the hyperbolic tangent and square-root functions mentioned in the footnote of section 5.2. For example, we can use the hyperbolic tangent formula, ν j = (ν f /2)[1 + tanh{γ(j − j max ′ )}], to construct the fraction of red molecules, R j = ν j − ν j−c , and then determine the location of its peak value by setting ∆R j = R j − R j−1 = 0. The result is j − j max ′ = (c + 1)/2, where j max ′ is the j-value for which ν j has its greatest slope. j max (without the prime) is where ∆ν j has its peak and differs from j max ′ by 1 unit according to j max = j max ′ +1. Consequently, j red −j max = (c+1)/2−1 = (c−1)/2, which is essentially the same as the starting value of 0.51(c − 1), that we found in Figure 19 for the RGB-profiles 2 .
Despite this agreement between the lag obtained for the RGB-model and the lag obtained for the hyperbolic tangent formula, there is an important difference. Namely, as c increases, the RGB-lag gradually deviates from its initially linear c-dependence due to the presence of the small quadratic term in Eq.(6.5). As shown in Figure 19 , the plot bends over and eventually saturates at 2 units of collision time when c becomes comparable to the width of the ν j profile. In contrast, the lag for the hyperbolic tangent model remains linear as c increases until the saturation occurs. This distinction reflects different curvatures of the RGB-profile above and below the point that ν j has its maximum slope, especially as c becomes large. We have already seen this asymmetry in the plots of ν j and ∆ν j in Figure 2 when the red molecules were permanently contagious, and in Figures 6 and 7 when c > 2. The hyperbolic tangent formula does not have this asymmetry and therefore gives a lag that remains linear as c increases until the saturation occurs.
We can understand why this asymmetry causes j red − j max to deviate from its initially linear dependence by referring to plots of R + G (black curve), G (green curve), and R (red curve) in Figure 15 , and considering how the black curve changes around j p , its point of maximum slope. For j ≥ j p , the slope can be approximated by s j = s p − a 2 (j − j p ), where a 2 is the curvature (i.e. the rate of change of the slope with respect to j) on the upper branch of the profile. Likewise, s j = s p − a 1 (j p − j) on the lower branch where j p ≥ j. We assume that both a 2 and a 1 are positive so that the slope decreases on both sides of j p , as it ought to do at a point of maximum slope.
The essential point here is that for each value of j, the red curve is produced by subtracting contributions from the upper branch of the black curve where j = j p + δ and the lower branch of the (identical, but shifted) green curve where j = j p − c + δ. Consequently, if c is small and the black and green curves lie close together, then the peak of the red curve will lie at the mid-point between the black and green curves where δ = c/2 and the slopes are both very close to the peak value s p . (As noted above, the condition ∆R j = 0 adds an extra half step, and the reference to j max , rather than j p , subtracts a full step, so that the result is j red − j max = (c − 1)/2, rather than c/2.) However, when c becomes large, the slope will be s p − a 2 (j − j p − δ) on the upper branch of the black curve and s p − a 1 (j p − j − c + δ) on the lower branch of the green curve. To find the value of j for which these slopes are equal, we simply equate those expressions and solve for j, obtaining j −j p = {a 1 /(a 1 +a 2 )}c. If the accelerations are equal, then j = j p +(c/2) (or j red − j max = (c − 1)/2 for the condition ∆R j = 0). But if a 2 > a 1 , then j − j p < c/2, and the location of the red curve falls behind the mid-point by the amount (c/2){(a 2 −a 1 )/(a 2 +a 1 )}. This asymmetry causes the lag to deviate from its initially linear dependence on c that we found in Figure 19 .
In summary, the 'herd immunity' occurs when red molecules are turning green as fast as they are being produced by the infection of blue molecules. At this time, the number of red molecules has reached its peak and the infection rate of blue molecules is slightly past its peak. In fact, it is convenient to think in terms of the time that the blue molecules are being removed most rapidly (and the red-plus-green molecules are being created most rapidly). This time occurs approximately half-way up the growth profile for red-plus-green molecules, which means that the red peak will occur somewhere in the time remaining between the middle and end of that profile. As Figure 19 showed, this time lag, j red − j max , depends on c, but it is relatively small and always less than 2 collision times, τ s . Thus, for serious social distancing and a value of c close to 1, the time, j max , will be greatly delayed, as will j red , which follows closely after it. However, the good news is that the immunity threshold, ν herd will be correspondingly small and so will the total number of infected molecules, ν f , when the epidemic is over.
Finally, it is important to recognize that the 'herd immunity' does not give the blue molecules any immunity against infection by the red ones. It is just a way of saying that the rate of infection has reached its maximum, and will be decreasing for the remainder of the ramp time. Although it may be reassuring to think of the green molecules as shielding blue molecules from red ones, it is not so reassuring to realize that red molecules provide the shielding when infected molecules are permanently contagious. Once the infection rate has reached its peak, the red molecules are so plentiful that some block others and prevent them from participating in the 'feeding frenzy' that will now continue at a diminishing rate. In the words of former baseball player, Yogi Berra (1998) , 'It ain't over 'til it's over' and all the red molecules are gone.
An interesting application of these calculations is to see what would happen if the amount of social distancing were relaxed suddenly before the virus has been completely eliminated. Figure 20 shows the result of suddenly changing c from c = 1.5 to 10 after 42 collisions when the fraction of infected molecules had increased to within 3.3% of its final value of 0.61, and the infection rate had decreased to about 14% of its peak height. This calculation was done for a population of N 0 = 10 5 , which means that only 0.033 × 0.61 × 10 5 ≈ 2000 more molecules would be infected if the social distancing were to continue at its current rate. What was the number of contagious molecules at that time? The fraction of contagious molecules is so small that we cannot obtain an accurate estimate from Figure 17 . And in that figure, the overlap between the curves with c = 1.25 and c = 1.5 make the estimate even more difficult. However, we can obtain a fairly accurate estimate from the enlarged plot of ν j in Figure 6 in combination with Eq.(5.1) for ∆ν j . In Figure 6 , the slope of the curve with c = 1.5 is 0.00657 per step (or collision) length when j = 40 − 42. For c = 1.5 (corresponding to 1.5 steps), this means that the change ν 42 − ν 40.5 = 1.5 × 0.00657 = 0.00986. Thus, for a total population of N 0 = 10 5 molecules, this corresponds to 986 contagious molecules after j = 42 collision times. This number is relatively small compared to the 0.6 × 10 5 molecules that had been infected at this time, as we can see in the bottom panel of Figure 17 . However, as we will see next, those 986 contagious molecules would soon do great damage.
As shown in the right panel of Figure 20 , at j = 42, the infection rate increased suddenly from ∆ν j ≈ 0.007 to ∆ν j = 0.13, which is 2.5 times the peak rate that occurred near j max = 32. By comparison, the fraction of infected molecules increased rapidly from ν j ≈ 0.59 (near their plateau value of ν j = 0.61) to ν j = 1. Thus, even when the infection rate had decreased to less than 1% (which is about 14% of its peak value), the sudden relaxation of social distancing caused a rapid burst of infection that effectively wiped out the 40,000 surviving blue molecules. Clearly, 986 contagious (red) molecules were sufficient to restart the infection at a very high rate (2.5 times the previous peak rate at j max = 32).
We have seen that c = 10 is equivalent to a permanently contagious virus with no social distancing. Now, we consider a presumably more realistic situation in which c changes suddenly from c = 1.25 to c = 2.50. Figure 21 shows the first case, in which the social distancing changes after 62 collisions when the fraction of infections is leveling off near ν = 0.35 and the infection rate is down to about 25% of its peak value. In this case, there is a surge in the infection rate, reaching approximately 0.04, which is about 2.0 times the peak height prior to the change. Also,the fraction of infected molecules changed from 0.35 to 0.80, corresponding to infection of 45,000 extra blue molecules in this population of N 0 = 1 × 10 5 . Of course, this increase in the number of infections is proportional to the size of the population and would be higher (or lower) if the population were higher (or lower). Figure 22 shows the second case, in which c also changes from c = 1.25 to c = 2.50, but after 48 collisions when the epidemic was near its peak. In this case, there is a much larger surge in the infection rate, reaching approximately 0.11, which is 5.5 times the peak height prior to reducing the amount of social distancing. Also, the fraction of infected molecules reaches about 0.89, corresponding to 54,000 more infections than would have occurred if the social distancing had been maintained at the value of c = 1.25. Thus, 9,000 more infections occurred when the change was at j = 48 near the peak of the infection rate than when the change was at j = 62 toward the end of the epidemic. As before, this enhancement was for a population of N 0 = 1 × 10 5 molecules, and would be proportionally more (or less) for a larger (or smaller) population. Making the change from c = 1.25 to c = 2.5 near the peak of the growth rate caused the final fraction of infected molecules to be ν f = 0.89. This value was larger than the value of ν f = 0.80 that was obtained in the previous example when the growth rate had fallen to 25% of its peak value. However, ν f = 0.89 is still less than the final value of ν f = 0.93 that would be obtained if the epidemic proceeded from the beginning with c = 2.5. Equivalently, ν f = 0.93 would have been obtained if the transition occurred for j ≤ 35 toward the start of the epidemic with c = 1.25. In fact, the more that the transition is delayed, the smaller the final fraction, ν f , becomes until the limit of ν f = 0.77 is reached for j ≈ 75. (Nevertheless, ν f = 0.77 is still larger than the value of ν f = 0.37 that would have been obtained if the change had not been made and the virus had progressed to completion with the original value of c = 1.25.) Likewise, the maximum growth rate decreases from ∆ν j = 0.16 to ∆ν j = 0.027 when the transition time moves from j = 35 to j = 75.
Now, let us see what happens when the value of c is decreased suddenly from c = 2.5 to c = 1.25, corresponding to an increase of social distancing. Figure 23 shows the result when the change occurs prior to the maximum in the infection rate while c = 2.5. Instead of continuing up to the final value ν f = 0.93, the fraction of infections levels off at ν f = 0.65, as shown in the left panel. Also, as shown in the right panel, the infection rate abruptly stops its rise and begins a more rapid decent toward 0. Figure 24 shows the effect of making the change at j = 22, after the peak in the infection rate. Again, the fraction of infections levels off quickly, but this time at ν f = 0.80, which is a larger value than we obtained by making the change earlier. Also, the infection rate falls suddenly, but too late to influence the height of the peak. So if we increase the amount of social distancing, we will obtain fewer total infections and a lower peak infection rate if we make the change early before the infection rate reaches its peak than if we wait until after the peak has occurred. This is opposite of what we found in the previous section where the amount of social distancing was suddenly decreased. In that case, it was necessary to wait until well after the peak to keep the final number of infections and the infection rate low. Thus, to minimize the total number of infections, social distancing should be applied as soon as possible when entering an epidemic and it should be removed as late as possible when coming out of the epidemic.
This paper describes a mathematical model for the spread of a virus through a population represented by the colliding molecules of a gas. In this model, uninfected molecules are colored blue, contagious molecules are red, and molecules that were infected, but have lost their ability to infect are colored green. The epidemic starts with a single red molecule entering a gas of N 0 − 1 blue molecules. The red molecule collides with a blue molecule, which becomes infected and turns red. Then the two red molecules collide with two more blue molecules making them red, and the process continues, rapidly increasing the number of infected molecules in the gas.
If infected molecules remained contagious forever, the fraction of infected molecules would continue to increase exponentially until about half of the molecules are infected. Then, with fewer blue molecules left, the infection rate would decrease and the epidemic would end when all of the remaining blue molecules are infected. No one has escaped. However, if an infected molecule loses its contagiousness after a finite time, τ c , then the spread of the infection is delayed and weakened, and the epidemic ends without infecting all of the blue molecules. The spreading rates, times, and magnitudes all depend on the ratio of the contagious lifetime, τ c , to the average time between collisions, τ s , the step rate used in the calculation. Thus, the key parameter is c = τ c /τ s .
In principle, we could fix τ s and let τ c vary, as if we were studying a variety of different diseases. However, I have regarded τ c to be fixed, and imagined that the collision time, τ s , varies with c due to different amounts of social distancing for a single disease.
After setting up the model with permanently contagious molecules, I performed the numerical computations and graphed the solutions, which included the evolution of the number of infected molecules as well as the infection rate. Also, for this case, it was possible to solve the difference equation for the spreading analytically, exactly reproducing the ramped transition between the initial and final states of infection as well as the peaked profile of the infection rate. For permanently contagious molecules, the final fraction of infections was ν f = 1, as expected. Also, the width of the infection-rate was about 4 units of τ s , the peak of the infection-rate profile was 0.25/τ s , and the location of the peak was at j ≈ log 2 N 0 +0.471, which is 17.08 for a population of 10 5 molecules.
In section 4, I dropped the assumption that the molecules were permanently contagious and derived the spreading equation for the fraction of infected molecules, ν j , for several values of the contagious lifetime parameter, c. Immediately, the plots showed that this more general solution reduced to the previous one if c was greater than about 10. So if the molecules can stay red for at least 10 collision times, they will wipe out the entire population of blue molecules, as if these red molecules remained contagious forever.
However, for smaller values of c, and especially values less than 2 and approaching 1, the results differed quantitatively from the previous solution for c > 10. The trends were the same, showing the same ramped transition from the initial state of uninfected blue molecules to the final state when the epidemic was over. However, the difference was that the ramp of rising infection had a lower slope and a longer duration than when the red molecules had an unlimited lifetime. Also, the ramp started at a later time than when the red molecules were permanently contagious. Most important, the decrease in slope was greater than the increase in duration, causing their product (and the final level of infection) to be less than 1. At the end of the epidemic, all of the red molecules were gone. The population consisted entirely of green molecules that had lost their contagiousness and blue molecules that had escaped infection. These changes became more pronounced as c fell below 2 and approached 1, corresponding to the so-called 'flattening of the curve' that happens with increased social distancing.
In the section on red, green, and blue molecules, we learned that the number of green molecules at a given time was equal to the total number of infected molecules at a time that is c steps earlier. This meant that the time history of green molecules is the same as that for all of the infected molecules, but delayed by c time steps. Both groups have the same threephase evolution with an initially flat distribution, followed by a steep ramp of rapidly rising infections, which levels off at the final fraction of infected molecules. The difference between these two groups gives the fraction of red molecules. When the shift, c, is relatively small, this difference is just c times the derivative of the profile of infected molecules. Consequently, for small c, the fraction of red molecules has a profile that is just c times the profile for the growth rate of all the infected molecules (expressed mathematically, ν j−1 − ν j−1−c = c∆ν j ). For c ≥ 3, this relation breaks down, and the rising phase of the profile for red molecules matches the ramped profile for all of the infected molecules. However, unlike the profile for the total number of infections (red plus green), the profile for red molecules reaches a peak and then returns to zero, leaving only green and uninfected blue molecules at the end.
Finally, in section 7, we examined the effect of changing the amount of social distancing during the epidemic. By suddenly removing all of the social distancing (i.e. increasing c to ∞, or, even to c = 10), all of the accomplishments disappeared and the fraction of infections rapidly shot up to the final value of ν f = 1. All of the the blue molecules became infected and the epidemic ended. For smaller changes either increases or decreases, the result was not quite so dramatic, but the trend was clear. To minimize the total number of infections, social distancing should be applied as soon as possible when entering the epidemic and it should be removed as late as possible when coming out of the epidemic.
The quantitative results, obtained from section 5 and appendices A-C, are:
1. The final fraction of infected molecules, ν f was well fit by the exponential expression The full-width at half maximum (FWHM) is somewhat less, ∼0.8w e , depending on the detailed shape of the profile.
3. The peak height of the infection rate was given by Here, σ c and ρ c can be obtained from Table 1 of Appendix A or from the best-fit approximations given in Eqs.(8.7a) and (8.7b) below. Substituting ν j = ν f , and setting j = j max − 1, we obtain
For c in the range (1,2) , σ c and ρ c are given by
where s = 1 + 4(c − 1). In addition, the best-fit approximations can be used over the full range of c:
For c > 3, ν f ≈ 1, σ c ≈ 1, and ρ c ≈ 2, so that Eq.(8.5) reduces to
which is the result when the red molecules are permanently contagious.
A point-slope approach gave an alternate expression for j max . In this approach, all of the linear plots of j max versus log 10 N 0 were assumed to pass through the same point, (log 10 N 0 , j max ) = (1.586, 6.153), so that the expression for j max became j max = 6.153 + 1 log 10 ρ c (log 10 N 0 − 1.586), (8.9) which reduces to j max − 0.884 = log N 0 log 2 = log 2 N 0 , (8.10) when c is large and ρ c ≈ 2.
6. The starting position or threshold value, j th , of the ramp of rapidly rising infection was defined by ν j th = 0.1ν f , where ν j th is obtained from Eq.(8.4) with j = j th , which is valid during the initial phase of the variation when ν j << 1. Consequently, j th becomes
Thus, 1/ log 10 ρ c gives the time interval between j th and j max − 1, the place that the curve of accumulated infections has its maximum slope. (By comparison, j max refers to the location of the maximum growth rate in the plot of ∆ν j versus j.)
7. The best-fit relation for the total fraction of infected molecules at the 'herd immunity' threshold was found to be ν herd = 1 − e −0.860(c−1) , (8.12) and the corresponding expression for the lag between the peak of the total infection rate, ∆ν j and the peak in the fraction of contagious (red) molecules, ν j − ν j−c , was found to be (8.13) until the lag saturated at c m = 8.81 and j red − j max = 2.
Now, we will use these calculations to understand the propagation of a virus. Remember that they are for an idealized model of a gas of molecules and not for a realistic population of people in our society. Also, this model and its calculations are only a first step in the experiment. A second step might include more complicated interactions including multiple collisions, clumps of molecules, and exchanges of molecules between separate populations. On the other hand, this simple model is general and ought to apply to any virus and population for specified choices of the parameters, c (equivalently, τ c ), collision time τ s , and N 0 .
We found that when the infected molecules retained their contagiousness indefinitely and there was no 'social distancing', the number of infections increased exponentially until about half of the molecules were infected. After that, the number of infections increased more slowly as the remaining molecules of the gas became infected. The infection rate, which we represented by ∆ν j , was essentially the derivative of the growth curve. It rose to its maximum value of 0.25N 0 /τ s in a time on the order of log 2 N 0 , which is the number of powers of 2 in the number N 0 (equivalently, 3.322 log 10 N 0 ). This means that the time to infect half the population would depend on the size of the population, and therefore would be longer for a metropolis than for a small community. In particular, for a small town of 10 4 molecules, it would take 13.3 collisions to infect half the population, which would be 2-4 weeks if the average time between collisions were 1 or 2 days. On the other hand, for a large city of 10 6 molecules, it would take 19.9 collisions and therefore about 3-6 weeks.
On the other hand, we found that the shape of the infection-rate curve was independent of N 0 , so that it had the same width and height for a small town and a large metropolis. In particular, its width (full width at half maximum) was about 4 collision times, which would be 4-8 days. Because 0.25N 0 /τ s × 4τ s = N 0 , this means that most of the molecules would be infected during this 4-step interval (of 4-8 days) around the peak of the distribution.
These results are a consequence of using Eq.(3.13a) to calculate ν j , the number of infected molecules expressed as a fraction of the population, N 0 . In this case, the final approach from ν = 0.5 to ν = 1 was slightly shorter than the initial rise from ν = ν 0 to ν = 0.5. Equivalently, the infection-rate curve was asymmetric around its peak with a steeper fall than rise, which means that the epidemic ends faster than it starts.
Next, I discarded the assumption that the infected molecules remain contagious indefinitely, and solved the problem for a range of contagious lifetimes, τ c , where τ c is c times the collision time (or step time), τ s . The result was a range of solutions that depend on the parameter, c. One way to think of this is to suppose that τ s is a constant and τ c varies, as if we were considering a variety of different diseases. Another approach is to suppose that τ c is constant, as it might be for a specific disease, and to consider what happens for a range of step times, τ s that would occur for different amounts of 'social distancing'. So, for the current pandemic, we can think of the c-dependence as an indication of the influence of 'social distancing' on the spread of the disease with c = 1 corresponding to perfect distancing and c = ∞ corresponding to the normal society with no social distancing.
The first result of these new calculations was to find that the plots of the number of infected molecules (expressed as a fraction of the total population, N 0 ) were similar to the original plot with no social distancing, showing a relatively long base level followed by a steep ramp to the final level where the infection stops. However, for increased amounts of social distancing, the curves took longer to reach the steep ramp, the slope of the ramp was lower, and the final level of infected molecules was smaller than in the absence of social distancing. Likewise the growth-rate profiles had smaller heights, larger widths, and were shifted to later times as the amount of social distancing increased. This corresponds to the popular term, 'flattening the curve'.
But how much flattening and delaying occurred? Let us assume that the contagious time is τ c = 14 days based on the time that an infected person is required to be quarantined for the current corona virus. As shown in Figure 7 , for c = 10 (corresponding to τ s = τ c /c ≈ 14/10 = 1.4 days between collisions), the curve peaks at j max = 17, corresponding to 1.4 × 17 = 23.8 days after the first contagious molecule entered the population of ∼10 5 molecules. For c = 1.25 (corresponding to 14/1.25 = 11.2 days between collisions), the curve peaks at j max = 48, corresponding to 11.2 × 48 = 537.6 days, or approximately 18 months. Referring to Figure 6 , we found that the fraction of infected molecules dropped from 1.0 to 0.38 (100,000 dropped to 38,000). So this social distancing saved 62,000 blue molecules, but required more than 1.5 years to do it. And this refers to the peak height, not the final infection rate, which reaches 0 around j max = 70, corresponding to 11.2 × 70 = 784 days, which is about 26 months. So the cost of those 62,000 blue molecules was more than 2 years of keeping apart and avoiding collisions for intervals of about 11 days.
The times are shorter for smaller populations and longer for larger populations, as we found in Eq. (5.10) and Figure 12 . For a small town of population N 0 = 10 4 , the peaks occur after times of 1.4 × 13.6 = 19 days for c = 10 and 11.2 × 35.5 = 397.6 days (about 13 months) for c = 1.25. For a metropolis with N 0 = 10 6 , the corresponding times are 1.4 × 20.2 = 28.3 days and 11.2 × 60.2 = 674 days (about 22 months). And these are only the times to reach the peak infection rates. The time required to bring those rates back to 0 is even larger. The number of blue molecules saved by social distancing is proportional to the populations, so for the smaller population of 10 4 , the number saved is 6,200, and for the metropolis of 10 6 , the number saved is 620,000.
As we have seen in Figure 6 , the growth of the numbers of infected molecules seems to go through three phases. The evolution begins with a long interval of low values, eventually changing to an upward ramp that ends in a plateau of final values. The width of the infection-rate profile is a rough measure of the duration of this upward ramp, and therefore an indication of how long the infection-rate remains at high values. As shown in Figure 7 and plotted in Figure 10 , the width, w e = ∆j, decreases as c increases. Therefore, the width increases as the social distancing increases. In fact, from Eq. (5.8a), one can show that the width, w e , is given approximately by w e = 5/(c − 1) for a large amount of social distancing with c 1.5. Thus, for the example with c = 1.25 (τ s ≈ 11.2 days), the width, w e , is 20 × 11.2 = 224 days, or about 7 months. Whereas, for our previous example with c = 10 (τ s = 1.4 days), the width was approximately 4 × 1.4 = 5.6 days. So when social distancing reduces the fraction of infected molecules from 1 to 0.38, it increases the duration of the active phase of the disease from about 1 week to 7 months. By social distancing, we trade a short duration of infection at a high rate for a long duration of infection at a low rate. This is the so-called 'flattening of the curve'. However, as Eq.(5.9a) shows, the benefit is to substantially reduce the final number of infected molecules.
It is also interesting to understand this process in terms of the numbers of contagious molecules (expressed as a fraction of the total population, N 0 ). For example, if a particular molecule managed to escape the virus and remain blue for a while, what would its chances of encountering a contagious molecule be at that time? Recall that the number of contagious molecules at a given step time, j, is equal to the number of infected molecules at that time minus the number of infected molecules c steps earlier at the time j − c. This difference depends on the size of the shift, c, relative to the width, w, of the infection-rate profile. For social distancing with c ≤ 2, the peak number of contagious molecules is approximately c times the maximum slope of the growth curve (which is the peak value of the infection-rate curve). So the resulting number of contagious molecules in the top panel of Figure 17 is just c times the corresponding peak heights in Figure 7 . Because c > 1 for an uncontained virus, this means that the peak heights for the numbers of contagious molecules are always greater than the corresponding peak heights for the infection-rate profiles. On the other hand, for a nominal amount of social distancing, the widths and locations of the peaks are the same for the contagious (red) molecules as they are for the infected (red plus green) molecules. Based on the numbers that we found in the examples above, this means that the surviving blue molecules must wait a few years before re-entering the pool if they wish to avoid collisions with contagious red molecules.
It is important to remember that this molecular model and its numerical calculations began as an attempt to determine the mathematical properties of the growth curve and growth-rate curve. The model was not intended to be a realistic description of the spread of the Covid-19 disease. Even the choice of 14 days for τ c was done arbitrarily to permit the calculation of some specific delay times that might occur in the evolution of the disease. The resulting delays of months and years might be different if another value of τ c were used in the calculation. Nevertheless, the reader should be pleased to gain some insight into how the process works, and not be too alarmed or complacent about the resulting numbers. Remember, these numbers refer to molecules, not people. Like many people, I was interested in how a virus might spread. Rather than looking for the answer in the published literature, I created the molecular model and did the calculations as if I were solving a puzzle in math or physics. Only after finishing these calculations, did I perform a Google search and find the three-component model of Kermack and McKendrick (1927) and a summary of that model by Weisstein (2004) . Of course, this was 'the tip of the iceberg', which led to more recent references, including those of Anderson and May (1979) , Jones and Sleeman (1983) , Smith and Moore (1996) , and others. The three components of the 'SIR-model' referred to individuals who were susceptible (S), infected and contagious (I), and recovered and not contagious (R) -the same components that I have called blue (B), red (R), and green (G), respectively. Therefore, we should probably call our molecular model the RGB model, after the RGB color coding used in electronic imaging and photography. Kermack and McKendrick described these populations with three differential equations, of which only two are independent because the total number of molecules is conserved. In our RGB terminology, R + B + G = N 0 , and the SIR equations are:
Here, 1/β corresponds to the collision time, τ s , and 1/γ is the (1/e)-lifetime of the red molecules in the absence of a source, as one can see by setting B = 0 in Eq.(9.1b). (Equivalently, 1/γ is the average lifetime of the red molecules.) In the RGB-model, the infected molecules remain contagious for a finite time τ c = cτ s before suddenly losing their ability to infect. Thus, c is equivalent to the average lifetime, 1/γ, in units of the collision time, 1/β.
As Weisstein pointed out, the key value determining the propagation or damping of these equations is the 'epidemiological parameter', R 0 = (β/γ)B, where B stands for the fraction of uninfected blue molecules, and β/γ is analogous to the quantity that we call c in the RGB-model. Thus, in Eq.(9.1b), dR/dt changes sign from positive to negative when (β/γ)B falls below 1, and the number of red molecules starts to decrease. Setting B = 1, we obtain R 0 = β/γ, which is the fraction of molecules that a single red molecule at the start of the epidemic will infect before it turns green. However, in the RGB-model, that red molecule would be responsible for 2 c − 1 infections because each infected molecule goes on to infect other blue molecules before the original red molecule suddenly turns green. With this distinction, c is analogous to R 0 , and the difference equations of the RGB-model are analogous to the differential equations of the SIR-model, as will be described in detail in Appendix D.
A. Solution For ν j << 1
We have seen that when c ≥ 10, the fraction of infected molecules, ν j , can be obtained from Eq. (3.7), which has an exact solution given by ν j = 1 − 2 −2 {j−(jmax−1)} , where j max is where the growth rate ∆ν j = ν j − ν j−1 = 1/4 and ν j−1 = 1/2. This is the solution when the infected molecules remain contagious for the duration of the epidemic, and all of the blue molecules eventually become red. Also, recall that when the non-linear term in Eq. (3.7) is neglected, the linearized equation becomes ν j = 2ν j−1 and its solution is ν j = 2 j ν 0 . This raises the question of what the initial-phase of the solution is for other values of c.
To find out, we linearize Eq. (4.3) by replacing the factor 1 − ν j−1 by 1 and writing
This linear equation has an exact solution composed of terms of the form, r j , where r refers to the roots of the equation, r c+1 − 2r c + 1 = 0, obtained by substituting r j into Eq. (A1). When c = 2, this equation is r 3 − 2r 2 + 1 = 0, whose three roots are r = 1, (1 ± √ 5)/2. Consequently, when c = 2, the exact solution to Eq. (A1) is
where the coefficients, −1, 1 ± 2/ √ 5 were determined by setting ν j /ν 0 ≈ 1, 2, 4 for j = 0, 1, 2 and solving the three equations simultaneously. The approximate numerical form of Eq.
which is quickly dominated by 1.894(1.618) j when j becomes larger than about 6. This is to be compared with 1.355(1.839) j when c = 3 and 2 j when c is greater than 10. Also, as we shall see in the next few paragraphs, the dominant terms are 3.732(1.366) j for c = 1.5 and 7.536(1.207) j for c = 1.25. These results are given in Table 1 for an expression of the form
where σ c is the coefficient of the dominant r j -term, and ρ c is the corresponding value of r.
As we would expect, this dominant term approaches 2 j as c becomes large. And as we will show next, it approaches {2/(c − 1)}c j as c → 1.
We are especially interested in the solution when c lies between 1 and 2. Recall that we used a linear interpolation for ν j−1−c in section 4. For c in the range (1,2), this linear interpolation is ν j−1−c = (2 − c)ν j−2 + (c − 1)ν j−3 . Consequently, Eq. (A1) becomes
We can use Eq.(A4) to compare this analytical solution of the linearized Eq.(A1) with the numerical solution of the non-linear Eq.(4.3). The results are plotted in Figure 25 for c = 10, 2, 1.5, 1.25 and N 0 = 1 × 10 5 . (When c = 10, the linearized equation is ν j = 2ν j−1 , whose solution is ν j = 2 j ν 0 .) Recall that the linearized equations involve the approximation 1 − ν j−1 << 1, which means that we expect the approximate solutions to be valid only for small ν j . In Figure 25 , we see that the solutions for c = 10 and c = 2 are fairly accurate for values of ν j less than about 0.2. However, for progressively smaller values of c, the dashed and solid curves begin to separate at smaller values of ν j , below 0.1 for c = 1.5 and then below 0.05 for c = 1.25. However, the values of ν f also decrease as c becomes closer to 1, so the discrepancies will be smaller when normalized to the values of ν f . 
Although the solution of the linearized version of Eq. (4.3) diverges from the solution of the complete non-linear equation, it remains accurate near the initial rising phase of the evolution where ν j << 1. Consequently, we can use the solution given by Eq.(A4) to estimate where the rapid increase of ν j begins. We begin by supposing that the starting point can be defined by a threshold value, arbitrarily taken to be ν j = 0.1ν f . The corresponding value of j, say j th , is determined by the equation
Next, we set j = j max − 1 in Eq.(5.11) and divide that equation by Eq. (B1) . This division eliminates σ c , ν f , and N 0 , and gives
from which we deduce that
Thus, we have an expression for the difference, (j max − 1) − j th , as a function of ρ c , which is given by Table 1 of Appendix A and the solution to the equation r c+1 − 2r c + 1 = 0 when c is an integer greater than 1. (Recall that j max − 1 and j max refer to different profiles. j max − 1 is the location of the maximum slope of ν j , which concerns us here, whereas j max is the location of the peak of ∆ν j .)
Note that our arbitrary threshold of ν th = 0.1ν f led to the shift (j max − 1) − j th = − log 10 0.1/ log 10 ρ c = 1/ log 10 ρ c because the logarithms were taken to base 10. But we could have used any base, including ρ c , in which case, the shift would have been − log ρc 0.1/ log ρc ρ c = log ρc 10. For permanently contagious red molecules, ρ c = 2 and the shift becomes log 2 10 = 3.322. For convenience, we could have chosen a threshold using a power of 2, like ν th = 2 −3 ν f , for example, so that the shift is just log 2 2 3 = 3 collision times. Thus, for permanently contagious red molecules, the threshold j th would be just 3 steps back from the point of greatest slope, which itself is given approximately by j max − 1 ∼ log 2 N 0 .
Also, note that because the shift is (j max −1)−j th = 1/ log 10 ρ c , the shift will double from 3 to 6 collision times when ρ c is changed from 2 to √ 2 = 1.414 (corresponding to c = 1.581 using Eq. (C1a)). And it will double again to 12 collision times on going to ρ c = 2 1/4 = 1.189 (c = 1.228), as one can verify in Figure 6 (but remembering that the threshold value is not just 0.1 (or 1/8 for the simpler case), but is 0.1ν f (or (1/8)ν f ). Likewise, for the nearly symmetric profiles, the end of the ramp occurs about the same number of collision times after the point of greatest slope.
C. Approximating the Linearized-Solution Parameters, σ c and ρ c Because the parameters, σ c and ρ c are important for understanding the initial behavior of ν j , I thought it would be useful to obtain approximate analytical expressions for those parameters, expressed as a function of c. Consequently, in Figures 26 and 27 , I have plotted the values of ρ c and σ c , respectively, obtained from Table 1 of Appendix A, supplemented by additional values chosen to fill in the range of c. In addition, I have superimposed dashed curves corresponding to the best-fit exponential solutions given in Eqs. (C1a) and (C1b). The general form of Eq.(C1a) was chosen so that it would give 1 when c = 1 and give 2 as c → ∞. That was much easier than choosing Eq.(C1b) so that it would reduce to 2/(c −1) as c → 1 and approach 1 as c → ∞. An additional exponential factor was necessary to weaken the 2/(c − 1)-term so that σ c would approach 1 more rapidly than 2/(c − 1) as c → ∞. Finally, note that σ c contributes to the initial behavior of ν j as a shift in the value of j. We can see this by writing
where j sh = log σ c / log ρ c whose logarithms can be taken to any base. In this case, Eq. (5.12) of the main text can be rewritten as
Thus, if c were greater than about 2, ν f and σ c would be approximately 1, ρ c would be approximately 2, and j sh would be negligible. Eq.(C3) would then become
which is the result for molecules that remain contagious indefinitely. Figure 28 shows values of j sh plotted as a function of c. For comparison, the dashed blue curve indicates j sh , calculated from Eq.(A9) and j sh = log σ c / log ρ c . Eq.(A9) is the approximate solution of Eq.(A5), which only applies when c lies in the interval (1,2). Nevertheless, the extension of this curve fits the data points for c > 2 reasonably well. A horizontal line of the form j sh = 0 would also be a good fit for c > 3. Consequently, we could also fit the data points using Eq.(A9) inside (1,2) and j sh = 0 for c > 3 outside. As I mentioned in the discussion section of this article, I constructed the RGB-model and solved its equations without any reference to the published literature. I just wanted to see if I could figure it out on my own, as if one of my scientific friends had presented it to me as a mathematics puzzle. It was only afterward, when I started looking through the published literature that I found the well-known SIR-model of Kermack and McKendrick (1927) .
Both models describe the spread of a virus through a large population of initially uninfected individuals, and track the population of the same three types of individuals, which I have called red (infected and contagious), blue (uninfected), and green (infected, but no longer contagious). However, the models take different approaches. Whereas the RGB-model is based on individual collisions and discrete equations, the SIR-model is continuous and is controlled by differential equations. The contrast is like comparing the discrete spreading by a random walk with the continuous spreading by an effective diffusion, reminiscent of the two ways of simulating the transport of magnetic flux on the Sun's surface (Leighton, 1964; Wang and Sheeley, 1994) .
Although both models suppose that the infected individuals do not remain contagious indefinitely, the RGB-model turns off the contagiousness abruptly after a specific time, τ c , whereas the SIR-model allows the contagiousness to decay exponentially at a rate γ. The models have two other relevant parameters, τ s , the step rate or average time between collisions of the molecules in the RGB-model, and β, the rate at which the contagious red molecules and uninfected blue molecules come together. In each case, the important quantity is the ratio of these respective parameters -τ c /τ s for the RGB-model and β/γ for the SIR-model. I have used c for τ c /τ s , and in this section, I will also use c for the ratio, β/γ.
We begin with the SIR equations Eqs.(9.1a), (9.1b), and (9.1c) in the discussion section of this paper. They are subject to the initial conditions: R(0) = ǫ = 1/N 0 , G(0) = 0, and B(0) = 1 − ǫ. This corresponds to our starting condition of one red molecule and N 0 − 1 blue molecules in the RGB-model. In each case, R + G + B = 1 in our normalized units.
Our first step will be to relate the quantities, R, R+G (which corresponds to ν j in the RGB-model), and B, to the variable, G. We do this by dividing Eq.(9.1a) by Eq.(9.1c) to obtain
Inserting c = β/γ and solving this differential equation, we obtain
where we have chosen the constant of integration to satisfy the initial conditions of the previous paragraph. And, of course, once R + G is known, B is also determined because their sum is 1:
It takes a little more algebra to obtain R. We begin by rewriting Eq.(D1) as
The solution to this differential equation is
where we have again determined the constant of integration so that the resulting equation satisfies the initial conditions of the first paragraph. Now, we have R, R + G, and B expressed in terms of G. Unfortunately, that is the end of our analytical solutions. Next, we use numerical integration to evaluate G as a function of time, t.
We begin by combining Eq.(9.1c) and Eq.(D5) as
from which we can obtain an implicit expression for G(t):
where again we have adjusted parameters, this time by setting γt = (γ/β)(βt) = (1/c)(t/τ s ) and expressed t in units of τ s .
As a practical point, we note that the presence of ǫ prevents the denominator of the integrand from vanishing at the lower limit of G = 0. It also provides the dependence on the population, N 0 = 1/ǫ. As a second point, we note that the upper limit, G(t), is bounded by a maximum value, which is 1 when c = ∞, but less than 1 for finite values of c. Because R = 0 at the end of the epidemic, it follows that the maximum value of G is also the maximum value of R + G, which we have called ν f in the previous sections. In this case, the maximum value of G is determined by setting the denominator of the integrand equal to zero.
This provides a method for calculating the value of ν f in this SIR-model. We just set 1 − G − (1 − ǫ)e −cG = 0, neglect ǫ compared to 1, and solve for c as a function of G. The result is
where G f is the final value of G and G + R, which corresponds to ν j in the RGB-model. Figure 29 shows this value of ν f (c) plotted as the black, dashed curve with the corresponding values obtained from our RGB-model, already shown in Figure 8 . The agreement is fairly good overall, and it is nearly perfect in the region where c < 1.4. As we shall see next, the agreement for small values of c is a general property of the two models. Equivalently, the disagreement increases as the exponential decay time for R increases in the SIR-model. Figure 8 ).
Next, we use Eq.(D7) to create a table of numerical values of G(t), and then use Eqs.(D2) and (D5) to create corresponding tables of R(t) + G(t) and R(t). The results are plotted in Figure 30 for values of c = 1.5, 2, 3 and 10. Referring to the corresponding RGB-plots in Figure 15 , we see again that the two models are similar for small values of c. In each model, for c = 1.5 and c = 2, the variation of G is essentially the same as the variation of R + G, but shifted c units later in time, and the variation of R is approximately c times the corresponding growth rate. This starts to change when c = 3, and by c = 10, R + G Figure 15 ). Note the scale changes between some of the images.
increases so fast that it reaches about 0.67 before G has begun its increase. By this time, R no longer clings to the R + G curve, and begins its descent like an exponential of the form e −t/τc while G shows the corresponding exponential increase, 1 − e −t/τc . For the RGB-model in Figure 15 , when c = 10, G continues to be a delayed version of R + G. Likewise, R rises with R + G until G starts to increase, and then falls to 0 as G increases to 1. Thus, when c = 10 or more, the ultimate decline of the number of red molecules reflects the steep rise of the number of green molecules in the RGB-model, but it reflects the more gradual, exponential rise of the number of green molecules in the SIR-model.
The exponential behavior of the SIR model is shown even more clearly for the plots with c = 20 and c = 50 in Figure 31 . In the SIR model, these plots show a sudden infection of the entire population followed by an exponential decrease of the red population as it gradually turns green on the time scale, τ c . Another approach is to display the SIR plots of R + G and d(R + G)/dt in the same way that we displayed ν j and ∆ν j in Figures 6 and 7. Figures 32 and 33 show those SIR plots. At first glance, the plots look remarkably similar. However, a closer look shows some differences. In Figure 32 , the asymptotic values of R + G are slightly less than the corresponding values of ν f , as we have already seen in the plots of ν f versus c in Figure 29 . Also, in Figure 32 , the 'knees' of the plots of R + G versus time are more rounded than those for the RGB-model in Figure 6 , especially for the larger values of c. This causes the SIR plots of the growth rate in Figure 33 to return to 0 more gradually than the corresponding RGB-plots in Figure 7 . Likewise the peak heights for the SIR growth rates are slightly less than those for the RGB-model, except when c < 2. Again, this shows the general tendency for the models to agree when c is small and there is a lot of social distancing. As a final part of our comparison between the plots in Figures 32 and 33 and in Figure 6 and 7, we note that the SIR plots of the growth rate reach their maximum values a few steps before the corresponding RGB growth rates reach their peaks. Figure 6 and show the delayed onset and the reduced final fraction of the infected population as c approaches 1.
In our discussion of 'herd immunity' in section 6.2, we defined the threshold to be the value of ν j where the number of red molecules, ν j − ν j−c , has its maximum value and is starting to decline. The resulting expressions are much simpler in the SIR-model than in the RGB-model. We begin with Eq.(9.1c) which tells us that dG/dt = γR, and therefore that d 2 G/dt 2 = γdR/dt. Thus, the herd-immunity condition, dR/dt = 0, implies that the slope, dG/dt, is a maximum when R is a maximum. (By way of contrast, in the RGB-model, R obtained its maximum value slightly before ν j−c obtained its maximum slope.) As Eq.(9.1b) indicates, dR/dt = βRB − γR, so that dR/dt = 0 also implies that B = 1/c and therefore -The value of ν herd , evaluated numerically for the RGB-model and from R + G = 1 − 1/c for the SIR-model. In each case, ν herd is the total number of infections (R + G) when the number of contagious red molecules reaches its peak. model and the RGB-model. The red curve is the root-mean-square best fit to the computed data points for the RGB-model using the exponential given by ν herd = 1 − e −0.860(c−1) (D10) and the dashed black curve is the SIR-expression for R + G = 1 − 1/c given by Eq. (D9b) . Although the two models show the same (normalized) number of infections for c < 1.5, the SIR-model shows a smaller number and a slower approach to 1 as c becomes larger.
In Figure 35 , we show the difference, j red − j max , between the times that R and R + G reach their respective peaks in the SIR-model. The j-values for each point were obtained -The time lag, j red − j max , between the location of the peak for infected molecules, R + G, and the peak for contagious molecules, R, expressed in units of the average time between collisions, τ s , and plotted as a function of the contagious parameter, c, for the SIR-model.
by evaluating the integral in Eq.(D7) for two different upper limits; for j red , we used the expression, ln(c)/c, given in Eq.(D9c); for j max , we used the root of the equation G+2e −cG = (c + 1)/c, which results from combining d 2 (R + G)/dt 2 ∼ G ′′ − c(G ′ ) 2 = 0 with the expression for G ′ = γ(1 − G − e −cG ) derived from Eqs.(9.1c) and (D5). The dashed red curve is the root-mean-square best fit to the data points, as given by the equation j red − j max = 0.946 = 0.806 ln(c + 0.221).
The data points used for this fit extend outward to c = 100, which has a j-value of j red − j max = 4.732. So although the lag does not increase rapidly with c, it nevertheless has not reached a limit by c = 100 (or even c = 1000, when the lag is 6.934). Like the corresponding lag determined for the RGB-model in Figure 19 , this best-fit logarithmic expression for the SIR-model has some curvature and increases very slowly with c. By comparison, 'purely antisymmetric' models, like the hyperbolic-tangent and square-root models described in the footnote to section (5.2) and in section (6.2), have lags that increase linearly as (c − 1)/2 before they saturate.
Finally, we note that the epidemiological parameter, R 0 , is often defined as the number of secondary infections that are produced in a homogeneous population of susceptible blue molecules by a single contagious red molecule before it loses its contagiousness and turns green. In the SIR-model, this works out to R 0 = β/γ, which we have called c. However, in the RBG-model, the number is 2 c − 1 because each infected molecule goes on to infect other blue molecules before the original red molecule suddenly loses its contagiousness. These secondary infections continue to infect blue molecules without the exponential weakening that gives the average decay time, 1/γ, in the SIR-model. As mentioned above, we can recover the classic SIR expression
by substituting c = R 0 in Eq.(D9b). However, the substitution of 2 c − 1 = R 0 in RGB Eq.(6.4) gives ν herd = 1 − e −0.860(c−1) = 1 − ( 2 1 + R 0 ) 1.24 .
These two equations give the value of ν herd (or R + G) when the number of red molecules reaches its peak. In Figure 34 we saw that these thresholds occur for different values of c. Now, comparing Eqs. (D12) and (D13), we can also see that the thresholds occur for different values of R 0 (except when c = R 0 = 1, so that ν herd = 0 and there is no epidemic).