# Ending the COVID-19 epidemic in the United Kingdom

L. C. G. Rogers
Statistical Laboratory, University of Cambridge 111 Statistical Laboratory,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WB, UK;
L.C.G.Rogers[AT]statslab.cam.ac.uk.

###### Abstract

Social distancing and lockdown are the two main non-pharmaceutical
interventions being used by the UK government to contain and control the
COVID-19 epidemic; these are being applied uniformly across the entire
country, even though the results of [3] show that the impact of the infection
increases sharply with age. This paper develops a variant of the workhorse
$SIR$ model for epidemics along the lines of [1],[7] [6], where the population
is classified into a number of age groups. This allows us to understand the
effects of age-dependent controls on the epidemic, and explore possible exit
strategies.

##  1 Introduction.

The global COVID-19 pandemic has swept through the nations of the world with a
frightening speed, and left governments struggling to cope with the situation.
The initial responses have been directed towards limiting the death toll and
ensuring that health services are not completely overwhelmed, as would be only
too possible with an infection that can grow by a factor of one thousand in a
month. As there is as yet no vaccine, no effective medication, and very
imperfect understanding of the parameters of the epidemic, efforts have been
directed towards containment, with decisions about return to normality being
left until later. Without vaccine or effective medical treatment, the only
remaining strategies would appear to be either a policy of contact tracing and
quarantining, or developing herd immunity. The first of these policies appears
to have been applied successfully in South Korea and Singapore, and is
generally regarded as the first line of public health defence. In the current
pandemic, most countries have quickly found themselves overwhelmed by the
scale and speed of the outbreak, and have been unable to apply contact tracing
as rigorously and universally as is needed for the method to work. When it
does work, contact tracing and quarantine will allow an outbreak to be snuffed
out before it spreads widely, but it will of course leave a large population
of susceptibles open to a new infection, so continuing vigilance is essential.
As we have seen contact tracing overwhelmed, the goal of this paper is to
explore the route to herd immunity, using age-dependent release from lockdown,
and a gradual relaxation of social distancing rules. In Section 2 we present
the model, which is in almost all respects a straightforward variant of the
standard $SIR$ epidemic model. The equations contain terms for the controls
which are available to modify the dynamics of the epidemic. The problem is a
control problem, and for this we have to define the objective, which we do in
Section 3. The issue is of course that we have a conflict between the obvious
cost of the numbers of citizens whose lives are ended prematurely, which is a
concern for the next few months; and the damage that an extended lockdown will
do to the economy, which will be a concern for many years if the aftermath of
the 2008 financial crash is any guide. In setting up the cost structure, some
relatively arbitrary (but hopefully reasonably realistic) assumptions have to
be made; these are not in any way essential to the approach, and can easily be
changed by any reader prepared to play with the Jupyter notebook posted
online222https://colab.research.google.com/drive/1tbB47uSGIA0WehY-
hvIYgdO0mpnZU5A8. Parameter values, or even the entire form of the costs, can
be changed by anyone with a little knowledge of Python. Experts in health
economics would doubtless be able to suggest values that better embody current
thinking, and before any of the results of this paper can be relied on, such
inputs will be necessary.

In Section 4 we briefly discuss the datasources used, and in Section 5 we
present the results of computation in various scenarios.

##  2 Model formulation.

A simple SIR epidemic model is too crude to allow us to model and control the
key features of the COVID-19 epidemic; many infected individuals are
asymptomatic, and the impact of the infection on different age groups is very
different. So we will break down the population into $J$ age groups, and let
$A_{j}(t)$, $I_{j}(t)$, $S_{j}(t)$ denote the numbers of $j$-individuals at
time $t$ who are (respectively) asymptomatic infected, symptomatic infected,
and susceptible. We will denote by $N_{j}(t)$ the total number of
$j$-individuals in the population at time $t$, and allow this to change
gradually with the influx of new births, visitors from other countries; this
is to model the possibility that new infecteds come in from outside and
reignite the epidemic.

The most basic form of the evolution is governed by the differential equations

$\displaystyle\dot{I_{j}}(t)$ $\displaystyle=$ $\displaystyle-\rho
I_{j}(t)+p\lambda_{j}(t)S_{j}(t)$ (1) $\displaystyle\dot{A_{j}}(t)$
$\displaystyle=$ $\displaystyle-\rho
A_{j}(t)+(1-p)\lambda_{j}(t)S_{j}(t)+\iota_{j}(t)$ (2)
$\displaystyle\dot{S_{j}}(t)$ $\displaystyle=$
$\displaystyle-\lambda_{j}(t)S_{j}(t)+\sigma_{j}(t)+\varepsilon\rho\\{\,I_{j}(t)+A_{j}(t)\,\\}$
(3) $\displaystyle\dot{N_{j}}(t)$ $\displaystyle=$
$\displaystyle\sigma_{j}(t)+\iota_{j}(t),$ (4)

where $\iota_{j}$ and $\sigma_{j}$ are known functions of time representing
the arrival of new asymptomatic infectives and susceptibles respectively333Any
symptomatic visitors attempting to join the population are assumed to be
prevented from entering.; and the final term on the right-hand side of (3)
allows for the possibility that removed infectives may not in fact be immune,
and some may return to the population ready for reinfection. The parameter
$p\in(0,1)$ appearing in (1), (2) is the probability that a susceptible
becoming infected is symptomatic; and the parameter $\rho>0$ is the recovery
rate. The infection rates $\lambda_{j}(t)$ are explicit non-linear functions
of the state of the system that will be discussed shortly, but, aside from the
terms involving $\lambda$, the evolution is linear. So if we stack the
variables into a single vector

$Z_{t}=[I(t),A(t),S(t),N(t)]$ (5)

the evolution (1)-(4) can be written as

$\dot{Z}(t)=MZ(t)+\Lambda(t)Z(t)+\eta(t),$ (6)

where $M$ is a $4J\times 4J$ constant matrix, $\Lambda$ is a simple matrix
whose entries involve the $\lambda_{j}$ in the appropriate entries, and $\eta$
is the vector of driving terms.

[Remark.The model (1)-(4) is the fluid limit of a Markov chain model in which
$\rho$ is the rate that an individual jumps from an infected state to the
removed state, and therefore the implicit (Markovian) assumption is that the
time spent in the infective state is exponentially distributed. This
assumption does not fit well with observation, so we can allow for different
distributions by the familiar trick of the method of stages (see, for example,
[2]), in which a infected individual passes through a number of exponentially-
distributed stages. In more detail, we can suppose that there are $K_{x}$
stages for the symptomatic infection, and that $I_{j}^{k}(t)$ is the number of
symptomatic $j$-individuals at stage $k$ of the infection at time $t$,
$j=1,\ldots,J$, $k=1,\ldots,K_{x}$. Making this change, the equation (1)
becomes the system

$\displaystyle\dot{I_{j}^{1}}(t)$ $\displaystyle=$ $\displaystyle-\rho
K_{x}I^{1}_{j}(t)+p\lambda_{j}(t)S_{j}(t)$ (7)
$\displaystyle\dot{I_{j}^{k}}(t)$ $\displaystyle=$ $\displaystyle-\rho
K_{x}(I^{k}_{j}(t)-I_{j}^{k-1}(t))\qquad(k=2,\ldots,K_{x})$ (8)

This corresponds to making the duration of symptomatic infection a sum of
$K_{x}$ independent exponentials each with mean $1/\rho K_{x}$, which has the
same mean as an exponential of rate $\rho$ but smaller variance. We could
similarly decompose the asymptomatic infections, and indeed by further
ramifications of the method of stages we could make the distribution of
infected time approximate any desired distribution. There is a good reason not
to take this too far, however; in the numerics, the differential equation has
to be solved many times. It is remarkable that this can be done in a
reasonable amount of time, but the more complicated the model, the slower this
step becomes and ultimately the computation will be too slow.

But however we do this, when we stack all the variables into a big state
vector $Z$, the evolution still has the form (6), and the appropriate form of
this is coded into the Jupyter notebooks.]

###  2.1 The specification of $\lambda$.

Each individual spends part of the waking day at home, and part of the waking
day outside444A more refined decomposition of time outside the home can be
made, but this does not really change the principles being explained here. .
We shall denote by $m^{O}_{ij}$ the mean number of contacts that an
$i$-individual has per day with $j$-individuals when outside the home; and by
$m^{H}_{ij}$ the mean number of contacts that an $i$-individual has per day
with $j$-individuals when inside the home.

It is important to understand that $m^{O}_{ij}$ is the mean number of contacts
that an $i$-individual has with $j$-individuals if everyone spends their
entire waking day outside the home. If the $i$-individual spends a fraction
$\varphi_{i}$ of the waking day outside the home, and $j$-individuals spend a
fraction $\varphi_{j}$ of the waking day outside the home, then the mean
number of contacts per day which an $i$-individual has with $j$-individuals
outside the home will be $\varphi_{i}m^{O}_{ij}\varphi_{j}$.

Each time an infected person has contact with someone, infection will be
transmitted with probability $\beta$, though of course this will only result
in a change if the person contacted was susceptible. Thus the overall rate at
which infection is passed in the outside world to $j$-individuals will be

$\lambda_{j}^{O}(t)S_{j}(t)=\beta\sum_{i}\\{\,A_{i}(t)+\delta
I_{i}(t)\,\\}\;\varphi_{i}(t)\,m^{O}_{ij}\,\varphi_{j}(t)\;S_{j}(t)/N_{j}(t),$
(9)

where $\varphi_{i}(t)$ is the fraction of time spent in the outside world by
$i$-individuals, and $\delta\in[0,1]$ is the proportion of symptomatic
infecteds who go into the outside world. In an ideal situation, this would be
close to zero, but many people with the infection get only mild symptoms and
may not self-isolate. The number of $j$-individuals at time $t$ is $N_{j}(t)$,
so the factor $S_{j}(t)/N_{j}(t)$ on the right-hand side of (9) is the
probability that a contacted $j$-individual is susceptible.

This may have the appearance of a conventional extension of an SIR model, but
one point to flag straight away is that the controls $\varphi$ enter
quadratically in the expression for the infectivity, whereas some authors use
only a linear dependence. This is erroneous.

What happens in the home is rather more difficult to deal with. We could
simply take (9) and change superscript $O$ to $H$, and $\varphi$ to
$1-\varphi$, but this would be incorrect, because an infected individual
outside may go through the day and infect a large number of people, but within
the home there are relatively few that could be infected, so the scope to
spread infection is much reduced - this is after all the rationale for locking
down populations.

Without the constraint on the number of infections imposed by the household
size, a single infected $i$-individual in the home would be firing infection
transmissions at $j$-individuals at rate

$\gamma_{i,j}(t)=\beta(1-\varphi_{i}(t))m^{H}_{ij}(1-\varphi_{j}(t)).$ (10)

Thus if $\tau$ is the mean infective time, during the period of infectivity
each infected $i$-individual in the home will fire a
Poisson($\gamma_{ij}(t)\tau$) number of infections towards $j$-individuals,
and therefore will fire a possible $Z\sim{\rm Poisson}(\gamma_{i}(t)\tau)$
number of infections towards all others, where
$\gamma_{i}(t)=\sum_{j}\gamma_{ij}(t)$. However, the number of infections that
can strike another individual cannot exceed $N-1$, where $N$ is the size of
the household in which the infected $i$-individual lives. Data from the Office
of National Statistics allow us to deduce the distribution555The data given
groups together all households with 6 or more members, so we assume that all
such households have exactly 6 members. of $N$. The mean number of individuals
at whom the infected $i$-individual fires infections is

$\mu_{i}(t)=E[Z\wedge(N-1)]=\sum_{k\geq 1}P(Z\geq k)P(N\geq k+1).$ (11)

This is the mean number of infections the infected $i$-individual could fire
at others during a period of mean length $\tau$, so we will simply suppose
that while infected

the $i$-individual in the home will be firing infections at rate
$\mu_{i}(t)/\tau$.

An infection fired at another will be supposed to strike a $j$-individual with
probability $p_{ij}^{H}(t)$ proportional to $m_{ij}^{H}(1-\varphi_{j}(t))$;
and given that it strikes a $j$-individual, the probability that a new
infection results will be $S_{j}(t)/N_{j}(t)$. Thus the analogue of (9) for
new infections of $j$-individuals in the home will be

$\lambda^{H}_{j}(t)S_{j}(t)=\sum_{i}\\{\,A_{i}(t)+\delta
I_{i}(t)\,\\}\;\frac{\mu_{i}(t)}{\tau}\;p_{ij}^{H}(t)\frac{S_{j}(t)}{N_{j}(t)}$
(12)

We combine these to give finally

$\lambda_{j}(t)=\lambda_{j}^{O}(t)+\lambda_{j}^{H}(t).$ (13)

[Remark. These assumptions represent a compromise; any honest treatment of
what goes on within households would appear to require a decomposition of the
population into groups according to different household compositions by age,
meaning the size of the statespace gets out of control - which would render
the calculation impractical. ]

It is worth emphasizing that there are just four controlling parameters in
this model: $\beta$, the probability that a contact results in a transmission;
$p$, the probability that an infected person is symptomatic; $\rho$, the
reciprocal of the mean infective time; $\varepsilon$, the probability that a
removed infective is still susceptible. Other values which are needed for the
calculations, such as the mean numbers $m^{O}_{ij}$, $m^{H}_{ij}$ of contacts,
can be found from published estimates.

##  3 Costs.

There are three components to the cost: the cost of lockdown, the cost of
social distancing, and the cost of deaths. We take them in turn.

###  3.1 Lockdown costs.

There will be a normal level $\bar{\varphi_{j}}$ for the proportion of time
spent by a $j$-individual outside the home; for the purposes of the
computations, the assumption here is that of the 112 waking hours of the week,
40 are spent in school or work, 20 are spent in social activities, and 52 are
spent at home, making $\bar{\varphi_{j}}$ equal to 60/112 for all age groups.

If a $j$-individual is locked down at level $\varphi_{j}(s)$ at time $s$, we
propose that the cost by time $t$ should be proportional to

$\int_{0}^{t}s\,\\{\,\bar{\varphi_{j}}-\varphi_{j}(s)\,\\}\;ds.$ (14)

For constant $\varphi_{j}$, this will be convex in $t$, which seems realistic;
a short lockdown (as for a public holiday) causes little damage, but as the
time away from regular work stretches on, the damage suffered increases more
rapidly, as businesses collapse and workers are made redundant. We will
consider strategies where for some $0<u<v$ (which may depend on $j$)

$\varphi_{j}(s)=\varphi_{j}(0)+\frac{\bar{\varphi_{j}}-\varphi_{j}(0)}{v-u}(s\wedge
v-u)^{+},$ (15)

where $\varphi_{j}(0)$ is the initial level of lockdown applied. At time $u$,
this starts to be relaxed in a linear fashion, being fully relaxed by time
$v$. Integrating (14) up to $v$ gives the cost of a $j$-individual being
locked down as

$C\,(\bar{\varphi_{j}}-\varphi_{j}(0))\,\\{\;u^{2}/2+u(v-u)/2+(v-u)^{2}/6\;\\}$
(16)

for some constant $C$. If we think that the social cost of an individual being
locked down for one year is $SC1$, then the constant of proportionality in
(16) is fixed so that the cost will be

$SC1\times\frac{\bar{\varphi_{j}}-\varphi_{j}(0)}{\bar{\varphi_{j}}}\times\\{\;u^{2}/2+u(v-u)/2+(v-u)^{2}/6\;\\}\;\frac{2}{365\times
365}.$ (17)

This then has to be summed over all the members of the population, with a
small reduction for retired people, who would presumably impact the economy
less if they were prevented from going out.

In the numerical implementation, we fix $v=u+5$; this reduces the number of
free parameters, and in any case reflects the realistic situation that once an
age group is freed from the lockdown restrictions they will quite quickly get
back to normal activity.

###  3.2 Social distancing costs.

Social distancing imposes costs; public transport will have to run at reduced
capacity, as will restaurants and theatres. But these costs are steady ongoing
frictions which do not keep people away from work for months on end. If the
social distancing policy means that at time $t$ the number of contacts outside
the home are reduced to a fraction $SD(t)\in(0,1)$ of the normal situation,
then we propose that the cost of this policy by time $t$ would be proportional
to666We make the same rule for all age groups, for simplicity; it might be
hard to maintain a social distancing policy that discriminated between age
groups.

$\int_{0}^{t}\\{\,1-SD(s)\,\\}\;ds.$ (18)

The form of $SD$ is available to choose, and in the computations we suppose
that $SD$ rises from the initial value $SD_{0}$ to final value 1 in a
piecewise-linear fashion through $NSD$ stages. This allows for the possibility
that social distancing could be gradually relaxed by opening more and more
classes of business or public assembly. Thus at some time $u_{0}$, $SD$ starts
to rise to the first staged value $SD_{1}$ at time $v_{0}$, where it remains
until $u_{1}$; from there it rises to the next staged value $SD_{2}$ at time
$v_{1}$, and so on. We suppose that the levels of the stages are equally
spaced, but this can easily be altered. If there is just one stage, the policy
starts at some value $SD_{0}$ and at time $u$ starts to rise linearly to 1 at
later time $v$, so the overall cost will be

$C(1-SD(0))(u+v)/2.$ (19)

By considering the effect of social distancing for a year, we fix the constant
to give cost

$SC1\times\theta_{SD}\times(1-SD(0))\times(u+v)/730,$ (20)

where $\theta_{SD}\in(0,1)$ expresses the pain of social distancing relative
to lockdown. In the calculations, we will allow the profile of social
distancing to be a more general piecewise-linear continuous functions,
permitting social distancing to be relaxed in stages and held at intermediate
values.

###  3.3 Death costs.

Making an estimate of the cost of the death of an individual is ethically and
procedurally quite a vexed issue. For the purposes of the calculations
reported in this paper and as default values used in the Jupyter notebook, the
assumption is that the cost of the death of an individual is proportional to
the expected number of further years that they would have lived; and that the
constant of proportionality is of the same order as $SC1$, the cost of an
individual being locked down for one year. So the code has a parameter
deathfactor which is used to scale $SC1$ for the calculations.

This is only part of the story however. We need to calculate the number of
deaths which will result from any particular policy, and this comes from the
calculated stream of removed symptomatic infectives, coming at rate $\rho
I_{j}(t)$ in age group $j$. Most of these will have recovered, but a
percentage of these will need hospitalization, and of those a percentage will
need critical care. The probabilities depend on the age of the patient, with
older patients at much higher risk; estimates are given in [8] and are quoted
in [3]. So we calculate the rate at which new critical care beds are required.
Based on an estimate for the number of days a critical care patient needs a
bed (taken to be 20 days), and knowing the total number of available critical
care beds, we can keep a running count of the number of critical care beds in
use, and then see how many of the incoming patients for critical care can be
accommodated. Those who can be accommodated survive with probability $p_{cc}$
(taken to be 0.5); those who cannot are assumed to die. It is assumed that
younger patients always take priority in allocating limited resources.

##  4 Data.

The code is built around the data assumptions in [3], who use nine age groups,
0-9, 10-19, 20-29, 30-39, 40-49, 50-59, 60-69, 70-79, and 80+. The
probabilities of hospitalization and critical care need for these age groups
are estimated by Verity et al. [8]. The population numbers for these age
groups come from the Statista web site
(https://www.statista.com/topics/755/uk/). The number of critical care beds in
England at the end of 2019 was around 4100, with around 11000 more planned at
the emergency Nightingale hospitals, so as an optimistic figure we took 12500
to be the number. The mean infectious period was taken to be 7 days, in line
with values in [3], but it seems this can be highly variable. Various values
were tried for $p$, the probability of an infected person being symptomatic,
but the baseline for this parameter was 0.3. Infectivity was taken to be 2.4,
in line with values proposed by [3], but again there appears to be quite a
wide range of possible values, as we see from [4]. The contact matrix values
$m^{O}$, $m^{H}$ are derived from [5]; they work with different age ranges, so
some pre-processing of their data had to be done; the code for this is
available from the author on request.

##  5 Computation.

The code for the calculations was written in Python, and is available in the
Jupyter notebooks for the reader to scrutinize and experiment with. The first
approach was to take the objective and minimize this using the Scipy routine
minimize, which acts as a wrapper to fourteen different methods, only a few of
which were possibilities due to the constrained nature of the problem. The
only routine which managed acceptable runtimes was SLSQP, but it turned out
that for virtually all randomly-chosen starting points, the end point was the
same as the start point; so this suggested the method which is used in the
Jupyter notebook, which is simply to randomly generate control rules of the
form discussed above, and focus on those which do best.

It is of course impossible to present more than just a few cases, but we can
explain what the default values for all the relevant parameters are, and then
show how the outputs vary as some of these get varied. As defaults, we have
taken

$\displaystyle{\tau}$ $\displaystyle=$ $\displaystyle 7$ $\displaystyle{\rm
Infectivity}$ $\displaystyle=$ $\displaystyle 2.4$ $\displaystyle p$
$\displaystyle=$ $\displaystyle 0.3$ $\displaystyle\delta$ $\displaystyle=$
$\displaystyle 0.05$ $\displaystyle\varepsilon$ $\displaystyle=$
$\displaystyle 0.05$ $\displaystyle{\tt phasein}$ $\displaystyle=$
$\displaystyle 5$ $\displaystyle K_{x}$ $\displaystyle=$ $\displaystyle
K_{y}=9$ $\displaystyle SD_{0}$ $\displaystyle=$ $\displaystyle 0.2$
$\displaystyle NSD$ $\displaystyle=$ $\displaystyle 4$ $\displaystyle SC1$
$\displaystyle=$ $\displaystyle 1e4$ $\displaystyle SD_{\rm end}$
$\displaystyle=$ $\displaystyle 0.9$ $\displaystyle{\tt CCstay}$
$\displaystyle=$ $\displaystyle 20$ $\displaystyle{\tt pcc}$ $\displaystyle=$
$\displaystyle 0.5$ $\displaystyle\theta_{SD}$ $\displaystyle=$ $\displaystyle
0.25$ $\displaystyle{\tt deathcosts}$ $\displaystyle=$
$\displaystyle[82.,72.22,62.5,52.86,43.33,34.,25.,16.67,10.]*SC1$

As initial values, we assume there are 50 asymptomatic infecteds in each of
the 9 age groups, and the initial vector $\varphi_{0}$ is

$\varphi_{0}=[5,5,10,10,10,10,5,5,5]*\bar{\varphi}/100$

The costs of lockdown are supposed to be less severe for the older age groups,
so we use

${\tt qcost}=[1,1,1,1,1,1,0.5,0.25,0.1]*SC1$

As mentioned before, we took the number Nbeds of critical care beds to be
12500. We ran the calculation for 900 days (except in the do-nothing example,
which ran for 200 days). We insisted that lockdown ends for all but the oldest
age group (80+) by day 400, and we imposed the condition that social
distancing reaches its end value $S_{\rm end}$ by day 840.

###  5.1 Base case: do nothing.

In this base case, we shall take $SD_{0}=0.995$ and
$\varphi_{0}=\bar{\varphi}$, which is the situation where no social distancing
and no lockdown happens. There are 86,000 deaths, and using the proposed cost
parameters, the cost of deaths is £16bn, the cost of social disruption is
£0.26bn. In this scenario, the epidemic is short and massive; as we see from
Figure 1, everything is over in about 100 days, with a peak number of new
daily cases for critical care of 25,000, and for hospital admissions of almost
120,000. Figure 2 shows that the critical care provision is completely
swamped, with nearly 70,000 critical care cases unable to get a critical care
bed and therefore dying without the necessary care. It is hard to imagine how
such a scenario could be thought acceptable.

Figure 1: Do nothing epidemic. Figure 2: Do nothing epidemic.

###  5.2 Lockdown and social distancing.

In this scenario, fairly tight lockdown and social distancing measure are
applied from the beginning and gradually relaxed. The costs of lockdown and
social distancing this time amount to £146bn, the death costs to £4.4bn, and
the total number of deaths was 16,100. The load of new cases is much more
manageable, with a peak of just over 2,200 new critical care cases, and about
16,000 new cases in all. All but the two oldest age groups are out of lockdown
within 100 days, but looking at Figure 4 we see that even after 800 days the
epidemic is far from over; once the oldest group is let out of lockdown and
social distancing has come to an end, the epidemic starts to take off again.
Most worrying here is that from 500 days on, every single critical care bed is
taken by a COVID-19 patient, and thousands of elderly patients needing a
critical care bed are unable to obtain one. This supports the proposition that
some form of social distancing will have to be maintained for a very long time
if no treatment or vaccine can be found.

Figure 3: Lockdown and social distancing. Figure 4: Lockdown and social
distancing.

###  5.3 Lockdown and social distancing with Infectivity = 2.8.

Next we see what happens if in fact the infectivity is higher than the middle
case value of 2.4 suggested in [3]. This time, lockdown and social distancing
costs remain at around £146bn, death costs are about £8.5bn, and the total
number of deaths is 38,700. The general picture looks like the previous
situation but more accentuated; there is a clear second surge after the oldest
age group is released from lockdown, and some 25,000 die without the critical
care they need as the hospitals are submerged with cases. This time,
saturation of the critical care facilities begins around day 600 and keeps
going. Even maintaining social distancing at 90% is not sufficient to hold
back the epidemic in the longer run.

Figure 5: Lockdown and social distancing, $R_{0}=2.8$ Figure 6: Lockdown and
social distancing, $R_{0}=2.8$.

###  5.4 Lockdown and social distancing with $p=0.2$.

If the probability that an infective is symptomatic is reduced to 0.2, the
outcome is improved, with death costs around £4.2bn, lockdown costs little
changed, and total deaths reduced to 17,500. Figure 7 shows two pronounced
peaks to the infection, the second again coinciding with the final relaxation
of restrictions. The critical care capacity only saturates at around day 650
this time. The epidemic is on a smaller and more manageable scale; peak
admissions to critical care are just over 600, peak hospital admissions just
over 4000. This is not surprising, since the proportion of those infected who
are symptomatic (and therefore open to possible complications) is lower.
However, there are more undetected asymptomatic infecteds going about in the
population, so the number of deaths is higher than in the base case; it is
clear from the pictures that towards the end the epidemic is beginning to get
out of control.

Figure 7: Lockdown and social distancing, $p=0.2$ Figure 8: Lockdown and
social distancing, $p=0.2$.

###  5.5 Lockdown and social distancing with $SD_{\rm end}=1.$

In this scenario, we find the costs of lockdown and social distancing to be
reduced to £124bn, death costs around £5.7bn. The number of deaths is 22,700.
What is most clear from Figure 10 is that from the time that the 70-79 age
group is released from lockdown around day 200, the epidemic gradually gets
more out of control, with critical care at full stretch from day 450 onwards,
and the numbers of older patients needing critical care and dying without it
growing 10,000 at the end of the run.

Figure 9: Lockdown and social distancing, $SD_{\rm end}=1.$ Figure 10:
Lockdown and social distancing, $SD_{\rm end}=1.$

Of course, it is only possible to display a few examples, which barely begins
to explore the diversity of behaviour that will arise as parameters are
varied. This is the purpose of the Jupyter notebook which can be found at

https://colab.research.google.com/drive/1tbB47uSGIA0WehY-hvIYgdO0mpnZU5A8

##  6 Conclusions.

This paper offers a simple model for the current COVID-19 epidemic; no account
is taken of spatial effects, which could make a big difference to any
conclusions. The treatment of the spread of the infection in the home is an
approximation, plausibly based perhaps, but still an approximation.
Nevertheless, the modelling assumptions are simple and compact, and permit
rapid exploration of possible responses of a non-pharmaceutical nature. The
calculations require assumptions about the initial state of the epidemic which
are essentially guessed. Even coming into the epidemic once it is under way,
it would be hard to get reliable values for the numbers of asymptomatic,
susceptible and immune people in the population, not least because there is at
the time of writing no test to determine whether someone has had the infection
and is now immune, and only a rather unreliable test whether an individual
currently has the infection. No account is taken of parameter uncertainty.
This is a natural area of enquiry, but at the moment it seems that the data
that would support strong conclusions is not yet available. As it seems that
the key parameters are known with very little precision, a highly detailed
model, or a sophisticated story about statistical inference may help less than
some rough exploration of possible parameter combinations; as the epidemic
evolves around the world, we will undoubtedly learn more of its
characteristics, which will allow us better to control it.

## Acknowledgements.

It is a pleasure to thank Josef Teichmann, Kalvis Jansons, Ronojoy Adhikari,
Rob Jack, Philip Ernst and Mike Cates for illuminating discussions. As
economists will insist on noting, they are not responsible for the errors
herein.

## References

  * [1] Roy M Anderson, B Anderson, and Robert M May.  Infectious diseases of humans: dynamics and control.  Oxford University Press, 1992.

  * [2] Andrew D Barbour.  Networks of queues and the method of stages.  Advances in Applied Probability, 8(3):584–591, 1976.
  * [3] Neil Ferguson, Daniel Laydon, Gemma Nedjati-Gilani, Natsuko Imai, Kylie Ainslie, Marc Baguelin, Sangeeta Bhatia, Adhiratha Boonyasiri, Zulma Cucunuba Perez, Gina Cuomo-Dannenburg, et al.  Report 9: Impact of non-pharmaceutical interventions (NPIs) to reduce covid19 mortality and healthcare demand.  2020\.
  * [4] Qun Li, Xuhua Guan, Peng Wu, Xiaoye Wang, Lei Zhou, Yeqing Tong, Ruiqi Ren, Kathy SM Leung, Eric HY Lau, Jessica Y Wong, et al.  Early transmission dynamics in Wuhan, China, of novel coronavirus–infected pneumonia.  New England Journal of Medicine, 2020.
  * [5] Kiesha Prem, Alex R Cook, and Mark Jit.  Projecting social contact matrices in 152 countries using contact surveys and demographic data.  PLoS Computational Biology, 13(9):e1005697, 2017.
  * [6] Rajesh Singh and R Adhikari.  Age-structured impact of social distancing on the covid-19 epidemic in India.  arXiv preprint arXiv:2003.12055, 2020.
  * [7] S Towers and Z Feng.  Social contact patterns and control strategies for influenza in the elderly.  Mathematical biosciences, 240(2):241–249, 2012.
  * [8] Robert Verity, Lucy C Okell, Ilaria Dorigatti, Peter Winskill, Charles Whittaker, Natsuko Imai, Gina Cuomo-Dannenburg, Hayley Thompson, Patrick GT Walker, Han Fu, et al.  Estimates of the severity of coronavirus disease 2019: a model-based analysis.  The Lancet Infectious Diseases, 2020.


11institutetext: Karlsruhe Institute of Technology, KaiserstraÃe 12, 76131
Karlsruhe, Germany 22institutetext: University of Hawai’i at MÄnoa, 2500
Campus Road, Honolulu, HI 96822-2217, USA 33institutetext: Alma Mater
Studiorum - UniversitÃ di Bologna, Via Zamboni, 33, 40126 Bologna BO, Italy
44institutetext: Istituto Nazionale Di Fisica Nucleare, Viale Berti Pichat,
6/2, 40127 Bologna BO, Italy 55institutetext: Oak Ridge National Laboratory, 1
Bethel Valley Rd, Oak Ridge, TN 37830, USA

# Modeling Distributed Computing Infrastructures for HEP
Applications00footnotetext: This manuscript has been authored in part by UT-
Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of
Energy (DOE). The publisher, by accepting the article for publication,
acknowledges that the U.S. Government retains a non-exclusive, paid up,
irrevocable, world-wide license to publish or reproduce the published form of
the manuscript, or allow others to do so, for U.S. Government purposes. The
DOE will provide public access to these results in accordance with the DOE
Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Maximilian Horzela 11<EMAIL_ADDRESS>Henri Casanova 22
<EMAIL_ADDRESS>Manuel Giffels 11<EMAIL_ADDRESS>Artur Gottmann 11
<EMAIL_ADDRESS>Robin Hofsaess 11<EMAIL_ADDRESS>GÃnter Quast
11<EMAIL_ADDRESS>Simone Rossi Tisbeni 3344<EMAIL_ADDRESS>Achim Streit 11<EMAIL_ADDRESS>FrÃdÃric Suter 55<EMAIL_ADDRESS>
###### Abstract

Predicting the performance of various infrastructure design options in complex
federated infrastructures with computing sites distributed over a wide area
network that support a plethora of users and workflows, such as the Worldwide
LHC Computing Grid (WLCG), is not trivial. Due to the complexity and size of
these infrastructures, it is not feasible to deploy experimental test-beds at
large scales merely for the purpose of comparing and evaluating alternate
designs. An alternative is to study the behaviours of these systems using
simulation. This approach has been used successfully in the past to identify
efficient and practical infrastructure designs for High Energy Physics (HEP).
A prominent example is the Monarc simulation framework, which was used to
study the initial structure of the WLCG. New simulation capabilities are
needed to simulate large-scale heterogeneous computing systems with complex
networks, data access and caching patterns. A modern tool to simulate HEP
workloads that execute on distributed computing infrastructures based on the
SimGrid and WRENCH simulation frameworks is outlined. Studies of its accuracy
and scalability are presented using HEP as a case-study. Hypothetical
adjustments to prevailing computing architectures in HEP are studied providing
insights into the dynamics of a part of the WLCG and candidates for
improvements.

##  1 Introduction

The physics programs in high energy physics (HEP) at the Large Hadron Collider
(LHC) would not be conceivable without the computing infrastructure provided
by the Worldwide LHC Computing Grid (WLCG) TDR-WLCG ; Update-WLCG . This
complex large-scale distributed computing system is crucial for the processing
and storage of the physics data produced by the collaborations at the LHC,
e.g., the Compact Muon Solenoid (CMS) CMS . With the start of the second phase
of the LHC program by the end of the decade, the so-called high-luminosity LHC
(HL-LHC), an increase by a factor of ten in the data production rate and a
proportional increase of the related computing demands are expected. Such a
dramatic increase of the storage and compute requirements will not be met
solely by traditional upgrades of the current WLCG resources. Hence, the
available resources will have to be utilized more efficiently. Identifying
changes in the design of an infrastructure as complex, heterogeneous and large
as the WLCG, so as to maximize performance and efficiency, is challenging.

To address the above challenge, we propose to study computing infrastructure
designs in simulation, which has proved useful in the past for determining the
original structure of the WLCG MONARC-Simulator ; MONARC-HLT . Simulation
enables the comparison of the execution of HEP applications on several
candidate compute and storage infrastructures without the need to build
testbeds and without disrupting the execution of workflows on WLCG production
systems. In this paper, we present a novel simulation tool and illustrate its
capacity to simulate HEP applications executions on complex large-scale
distributed computing systems accurately.

##  2 Simulation Tool

We have developed a simulator in C++ to enable the simulation and performance
analysis of HEP applications using two state-of-the-art discrete-event open-
source simulation frameworks: SimGrid Casanova2014 ; simgrid-web  and WRENCH
wrench ; wrench-web . The simulation models in both frameworks have been
validated 10.1145/2517448  and make it possible to achieve sensible trade-offs
between simulation accuracy and simulation scalability (i.e., speed and memory
footprint). Hereafter we describe these frameworks, the simulator’s
implementation, the simulators calibration and validation procedure, and the
computational complexity of the simulator.

###  2.1 Supporting Simulation Frameworks

SimGrid provides low-level simulation abstractions, implemented in C++, for
the development of simulators of distributed computing systems. Some of these
abstractions are used to describe hardware platforms, which comprise network,
compute, and storage resources. Network resources consist of network links,
defined by latencies and bandwidths. Compute resources consist of multicore
processing units, defined by numbers of cores and core compute speeds. Storage
resources consist of disks, described by read and write bandwidths. Routes can
be defined between these resources, as vectors of network links, to define the
platform’s network topology. Other abstractions are then used to describe
simulated processes, or actors that execute on this platform. Each actor can
execute arbitrary C++ code and spawn activities (i.e., communication,
computation, or I/O) that utilize hardware resources. A key feature of SimGrid
is that it formulates the problem of computing the rate of progress of all
simulated activities at a given time as a unified constrained optimization
problem. The constraints involve the activities’ rates of progress, scaling
factors, and hardware resource capacities. The objective function corresponds
to Max-Min fairness, i.e., maximizing the minimum rate of progress over all
activities, where these rates of progress are considered to be real numbers.
This problem must be solved each time an activity starts or completes and
SimGrid implements a highly optimized solver that is invoked repeatedly
throughout the simulation’s execution. This optimization problem formulation
makes it possible to accurately model complex real-world behaviors with low
computational complexity 10.1145/2517448 .

WRENCH builds on SimGrid to ease simulator implementation wrench . It comes
with already implemented abstractions of commonplace components of such
systems, which can be used as building blocks for implementing simulators of
complex systems, while benefiting from the accurate and scalable simulation
models provided by SimGrid. WRENCH provides implementations for several kinds
of services such as compute services (e.g., bare-metal servers, batch-
scheduled clusters, HTCondor deployments, cloud infrastructures) and storage
services (e.g., simple file servers, proxy file services with caching
capabilities, XRootD deployments). The user defines a hardware platform and
defines which services are deployed on that platform and on which resources.
The user also implements one or more simulation controllers, i.e., processes
that interact with the deployed services to execute some application workload
of interest. This interaction is implemented through each service’s API. For
instance, WRENCH provides the notion of a job that comprises tasks with
arbitrary dependencies, can be submitted to a compute service, and eventually
succeeds or fails.

###  2.2 Simulator Implementation

At the onset of this project, WRENCH provided all abstractions necessary to
simulate the execution of batch computation, where each job simply consists of
reading input files, performing computation, and writing output files.
However, such batch jobs make up only a small fraction of the workloads in
HEP. Most jobs, instead, are streaming jobs that pipeline input reading,
computation, and output writing with a granularity defined, e.g., by the
XRootD block size. To support the simulation of streaming jobs, we implemented
in WRENCH the concept of a custom action that allows the user to define a
job’s behavior, such as streaming, based on an arbitrary C++ function.

Besides streaming jobs, other key behaviors of real-world executions that we
wish to simulate are data locality and caching. Specifically, when attempting
to read input data, a job checks a list of cache locations where the data may
be available. If found, then the data is obtained from that location.
Otherwise, the data is obtained from authoritative remote storage locations
that hold input data at the onset of the execution. In this case, the data is
then copied to the caches, which will be able to answer future requests for
access to this same data. When a cache’s storage capacity would be exceeded,
data is evicted based on some classical eviction policy, such as Least
Recently Used (LRU). We implemented this behavior by using another custom
action. Note that this caching feature is useful in general and, as such, has
been implemented in WRENCH as part of the storage service.

WRENCH provides features for monitoring of the simulation. To avoid a too
large amount of monitoring data from large-scale executions of our simulator,
we implemented our own lightweight simulation monitoring component to track
relevant information for our purpose: job identifiers, job execution times,
data cache hit rates, and data transfer times.

As explained in the next section, we simulate the execution of multiple
concurrent workloads. A distinct controller is in charge of the execution of
each of these workloads. Each controller interacts with storage services
deployed on the platform including XRootD, and with an HTCondor compute
service that provides access to and schedules jobs on a set of bare-metal
compute services deployed on the platform’s compute resources.

###  2.3 Simulator Usage

The simulator takes as input: (i) a description of the compute platform; (ii)
a specification of the workloads that execute on this platform; and (iii) a
specification of the data present in caches at the onset of the execution. We
configure a workload as a set of individual jobs, each executing pipelined
input reading and computation based on input file sizes and a data transfer
block size, and transferring an output file to a storage. Each workload is
characterized by various parameters, such as the number and type of jobs, the
submission time, probability distributions for the number and size of input
files, the required CPU cores and memory, a computational volume per byte of
input data, and the size of the output file.

The above input fully characterizes the execution to be simulated, including
the scheduling of the jobs and the execution of their underlying activities,
which trigger computations, I/O operations, and network communications.
Furthermore, input files are cached when a cache storage location is available
in a predefined vicinity of the machine that executes a job, and LRU files are
evicted from that cache when full. When a job completes, its output files are
written to the predefined storage locations. When all jobs have completed, the
simulation terminates, and all relevant information to study each individual
job is provided.

###  2.4 Calibration & Validation

The simulation models implemented in the simulator are meant to mimic the
empirically observed behavior of application executions on real-world
platforms. So as to capture a range of possible execution scenarios, these
models come with free parameters that have to be tuned to best match real-
world behaviors. Picking appropriate values for, or calibrating, these
parameters was performed based on ground-truth data. This data comprises a set
of collected execution traces of jobs executions on a controlled, dedicated,
and distributed test platform, where the jobs and the platforms have known
characteristics. The calibration process consists in quantitatively comparing
these real-world traces to traces generated by the simulator, and adjusting
parameter values to achieve best possible match.

Once the simulator is calibrated it is validated against ground-truth data
other than that used for performing the calibration. When the simulated
execution traces match those in the ground-truth data, we conclude that the
simulation models, and their calibration, are valid. A key performance metric,
that is computed from real-world or simulated execution traces, is the job
execution time, on each machine in the platform. This metric depends highly on
the fraction of input files that are initially in cache at the beginning of
the execution (which we call the hitrate). Fig. 1 shows job execution times
vs. hitrate, in the real-world and in simulation for a platform with three
worker nodes. All worker nodes have the same architecture but two of them have
12 cores and the third has 24 cores.

\begin{overpic}[width=212.47617pt]{figures/Computing/SimGrid/Validation/Measurements/hitrateWalltime_validationSGbatch_GW1G_HDDCache-1.pdf}
\put(50.0,48.0){\includegraphics[width=39.0242pt]{figures/Computing/SimGrid/Validation/SGBatch_validation_platform.pdf}}
\end{overpic}\begin{overpic}[width=212.47617pt]{figures/Computing/SimGrid/Validation/Simulation/hitrateWalltime_1G_HDDCache_XRD100M_BInf-1.pdf}
\put(50.0,54.0){\includegraphics[width=21.68231pt]{figures/Computing/SimGrid/simgrid_logo.pdf}}
\put(60.0,50.0){\includegraphics[width=21.68231pt]{figures/Computing/SimGrid/wrench-
logo-vertical.png}} \put(51.0,48.0){\small{DCSim}} \end{overpic}

Figure 1: Job execution times per execution machine (worker node) vs. the
hitrate in the real world (left) and the simulation (right). Plots show
medians and 2.5- and 97.5-percentiles of all executed jobs on each execution
machine.

The calibrated simulation is able to reproduce the observed medians of the job
execution times in the validation data without any adjustments to the
simulation parameters. However, as hitrate values increase, the simulated
spread in the execution times is increasingly underestimated. This is expected
since the simulation does not capture all aspects of the real-world system. In
particular, the complex I/O access patterns for the hard disk drives that
store cache data are not simulated. These I/O accesses, due to seek times and
random accesses, are non-deterministic. The simulator, however, models I/O
access times simply based on latencies and read/write bandwidth rates of the
storage resources, which leads to a much lower variance of job execution times
than in the real world. For flash drive storage the seek times would be both
smaller and more deterministic, which could be simulated with higher accuracy.

###  2.5 Computational Complexity

To assess the space- and time-complexity of the simulator the maximum memory
consumption and execution time of the simulation as functions of the size of
the simulated platform are studied. Fig. 2 shows the scaling of the
simulation’s memory consumption and execution time.

\begin{overpic}[width=212.47617pt]{figures/Computing/SimGrid/Scaling/hitrateWalltime_1G_HDD_XRD100M_BInf_scaled-1.pdf}
\put(50.0,54.0){\includegraphics[width=21.68231pt]{figures/Computing/SimGrid/simgrid_logo.pdf}}
\put(60.0,50.0){\includegraphics[width=21.68231pt]{figures/Computing/SimGrid/wrench-
logo-vertical.png}} \put(51.0,48.0){\small{DCSim}} \end{overpic}

Figure 2: Simulation execution time and maximum memory consumption vs. the
maximum number of active job slots, a proxy for the platform size (left); and
simulated job execution times vs. hitrate with mitigated time complexity
(right).

The simulation execution time is driven by the effective size of the platform,
that is the maximum number of jobs that can execute concurrently in active
_job slots_. A superlinear increase in execution time as platform size
increases is observed.

Given the above, a mitigation procedure is needed to enable the simulation of
large platforms under a time budget. By scaling the number of cores on each
simulated machine by a constant factor the effective number of active job
slots in the simulation of a large platform, and therefore the simulation
execution time, can be decreased superlinearly. To conserve the dynamics of
the simulation the relative proportions of all other platform parameters are
scaled by the same factor. The mitigation procedure is applied on the
simulation of the validation platform presented in Section 2.4. Comparing the
simulation traces of the full simulation, visualized in Fig. 1, and mitigated
simulation, visualized in Fig. 2, no significant deviation can be observed.
Consequently, the mitigation procedure makes the simulation scalable while
conserving the validity of the original simulation.

##  3 Proof-of-Concept: Diskless Tier 2

To showcase the applicability of the presented simulator for modeling the
performance of such architectures, we study two aspects of a future candidate
design for a subset of the WLCG infrastructure proposed by the German HEP-
computing community in KET-Pos .

###  3.1 Simulated Workload and Platform

The simulated workload characteristics for these studies are obtained from job
monitoring data gathered by the CMS collaboration of jobs executed at the WLCG
Tier 1 site at KIT and Tier 2 site at DESY in a time frame from February 24th
to March 7th 2023. From this data, and for each job, we extracted information
about computational volume, amount of read and written data, consumed memory,
and number of CPU cores used.

The jobs are classified as Analysis, ReadoutSim, Processing, Merge, and Other
jobs, which make up approximately 49.7%, 39.5%, 7.6%, 3% and 0.2% of the
workload, respectively. Analysis jobs contain all jobs submitted by individual
users and therefore shows a large spread in the distributions of their job
characteristics. ReadoutSim jobs process the simulation of the readout
electronics of the CMS detector CMS  and the digitization of the simulated
signals. Processing jobs process the reconstruction from the raw digitized
detector signals to analysis objects usable in analyses. Merge jobs combine
multiple data files created in multiple jobs into single files simplifying
subsequent data processing. Other jobs contain all generation, simulation,
reconstruction, cleanup, and logging executions that did not occur often
enough in the monitored time frame to result in a significant class.

From the individual characteristics of the jobs in each class we derived a
one-dimensional binned probability distribution for each characteristic
quantity. Then, we sampled the characteristics of the simulated workload from
these probability distributions while keeping the proportions relative to
those observed in the monitoring data. We neglected the correlations between
the probability distributions of the individual characteristics.

The simulated platform consists of two hypothetical sites approximately
representing the real-world WLCG Tier 1 and Tier 2 sites at KIT and DESY. Each
site is modeled as its own local network of worker nodes and storage
resources, and the two sites are interconnected via WAN. The exact platform
parameters are listed in Table 1.

Table 1: Characteristics of the platform input into the simulation for the PoC study. Characteristic | Tier 1 | Tier 2'
---|---|---
Compute | 42,000 cores | 20,000 cores
Storage Bandwidth | $80\text{\,}\mathrm{Gbit}\text{\,}{\mathrm{s}}^{-1}$ | n/a
Cache Bandwidth | n /a | $80\text{\,}\mathrm{Gbit}\text{\,}{\mathrm{s}}^{-1}$
LAN Bandwidth | $2\times$100\text{\,}\mathrm{Gbit}\text{\,}{\mathrm{s}}^{-1}$$ | $40\text{\,}\mathrm{Gbit}\text{\,}{\mathrm{s}}^{-1}$
WAN Bandwidth | $100\text{\,}\mathrm{Gbit}\text{\,}{\mathrm{s}}^{-1}$

The two sites are considered to be isolated from the greater context of the
WLCG, meaning that there is no additional load on the WAN caused by workloads
running on other sites.

###  3.2 Study 1: Grid Storage Replacement

The first simulation study evaluates the performance of a WLCG Tier 2 site
that is not a full-fledged storage site (which would require maintenance and a
personnel budget) but instead serves only as a simple stand-alone cache with
no availability or persistence guarantees (Tier 2'). All input files in the
simulation are served by the grid storage at the Tier 1 site. A defined
fraction of input files are additionally staged at the cache at Tier 2' at
start of the simulation. This makes it possible to evaluate the performance of
the workload execution dependent on the fraction of cached files (hitrate)
with repeated simulations for multiple hitrate value between zero and one. A
hitrate $=1$ corresponds to the nominal case of a Tier 2 site following the
WLCG definition. A hitrate $<1$ corresponds to the hypothetical scenario where
the managed storage at the Tier 2 is reduced to a cache with a performance of
this cache parameterized via the hitrate. At a hitrate $=0$, there is no cache
and all input files are served by the Tier 1 storage.

To estimate the performance of each site we measure the distribution of the
fraction of time spent processing over the total job execution time (CPU
efficiency) of the executed jobs. The obtained distributions in simulation are
visualized in Fig. 3(a).

(a) Study 1

(b) Study 2

Figure 3: Job CPU efficiency vs. the hitrate for the scenarios in Section 3.2
(left) and Section 3.3 (right). Results shown as standard box plots with the
median, 25- and 75-quartiles, minimal and maximal single values within the
interval spanned by the lower (upper) quartile minus (plus) 1.5 times the
interquartile range, and outliers.

It can be observed in Fig. 3(a) that a replacement of the grid storage with a
cache leads to a decrease in the performance of both sites. For hitrates $<1$
the bandwidth to the grid storage at the Tier 1 is shared by the jobs on both
sites which increases the idle time of the CPUs due to the jobs waiting for
data to process. For hitrates $\gtrapprox 0.5$ the performance of the Tier 2'
stabilizes due to the internal network throttling the execution. The Tier 2
load on the Tier 1 storage still limits the performance of the Tier 1. To
overcome the performance drop, an increase of the Tier 1 storage bandwidth
needs to be considered.

###  3.3 Study 2: Cache Addition at High Latency Site

In this study we explore the impact of replacing a dedicated Tier 2 site by a
hypothetical supercomputer. Hence, the study presented in Section 3.2 is
repeated but the WAN bandwidth between the sites is reduced by a factor of
ten, which corresponds to a scenario in which the supercomputer does not have
a high WAN bandwidth. This drastic reduction in the WAN bandwidth changes the
dependence of the performance of both sites on the hitrate drastically.

The effect of the lower network bandwidth between the sites is that the
execution time of jobs at Tier 2' is dominated by the time spent waiting for
input data for hitrate $\lessapprox 0.7$. Above this value, job executions at
Tier 2' are once again throttled by the relatively low internal network
bandwidth at that site. This leads to less load introduced by these jobs on
the grid storage at Tier 1. As a consequence, job executions at Tier 1 are
only marginally affected, but the performance of job executions on Tier 2' are
significantly improved by high hitrates.

##  4 Conclusions

We have proposed using simulation to predict the performance of current and
hypothetical architectural designs for distributed computing platforms such as
the WLCG. To this end we have implemented, calibrated, and validated an open-
source simulator Horzela2023  based on the SimGrid and WRENCH simulation
frameworks. The simulator’s flexible configuration and its use of scalable
simulation models make it possible to simulate and analyze the execution of a
wide range of application workloads on large-scale computing platforms under
an acceptable time budget. These capabilities were demonstrated in two case
studies, which involve systems with collectively more than 60,000 cores and
$10\text{\,}\mathrm{PB}$ of storage integrated into a distributed network, and
applications transferring and processing petabytes of HEP data.

A future direction is the simulation of larger distributed computing systems
up to the size of the WLCG and beyond. This will require reducing the
computational complexity of the underlying simulation models while preserving
their accuracy.

## Acknowledgements

This work was supported by the German Federal Ministry of Education and
Research (project FIDIUM 05H21VKRC2) and the Institute of Experimental
Particle Physics, the Steinbuch Centre for Computing and GridKa at the
Karlsruhe Institute of Technology, Germany; by National Science Foundation
awards 2103489 and 2106059; and by Laboratory Directed Research and
Development Strategic Hire funding No. 11134 from Oak Ridge National
Laboratory, provided by the Director, Office of Science, of the U.S.
Department of Energy.

## References

  * (1) K. Bos, N. Brook, D. Duellmann, C. Eck, I. Fisk, D. Foster, B. Gibbard, C. Grandi, F. Grey, J. Harvey et al., CERN-LHCC-2005-024 (2005)
  * (2) I. Bird, P. Buncic, F. Carminati, M. Cattaneo, P. Clarke, I. Fisk, M. Girone, J. Harvey, B. Kersevan, P. Mato et al., CERN-LHCC-2014-014 (2014)
  * (3) The CMS Collaboration, Journal of Instrumentation 3, S08004 (2008)
  * (4) I.C. Legrand, H.B. Newman, _The MONARC Toolset for Simulating Large Network-Distributed Processing Systems_ , in _Proceedings of the 32nd Conference on Winter Simulation_ (San Diego, CA, USA, 2000), p. 1794â1801, ISBN 0780365828
  * (5) H.B. Newman, I.C. Legrand, AIP Conference Proceedings 583, 164 (2001)
  * (6) H. Casanova, A. Giersch, A. Legrand, M. Quinson, F. Suter, Journal of Parallel and Distributed Computing 74, 2899 (2014)
  * (7) SimGrid Team, _SimGrid: Versatile Simulation of Distributed Computer Systems_ , [Online, accessed on 08.28.23] https://simgrid.org/ (2022)
  * (8) H. Casanova, R. Ferreira da Silva, R. Tanaka, S. Pandey, G. Jethwani, W. Koch, S. Albrecht, J. Oeth, F. Suter, Future Generation Computer Systems 112, 162 (2020)
  * (9) WRENCH Team, _WRENCH: Cyberinfrastructure Simulation Workbench_ , [Online, accessed on 08.21.23] https://wrench-project.org/ (2023)
  * (10) P. Velho, L.M. Schnorr, H. Casanova, A. Legrand, ACM Transactions on Modeling and Computer Simulation (TOMACS) 23 (2013)
  * (11) Komitee fÃr Elementarteilchenphysik, _Perspektivpapier der Teilchenphysiker:innen in Deutschland_ , [Online, accessed on 08.28.23] https://www.ketweb.de/sites/site_ketweb/content/e199639/e312771/KET-Computing-Strategie-HL-LHC-final.pdf (2022)
  * (12) M. Horzela, H. Casanova, A. Gottmann, S. Rossi Tisbeni, _HEPCompSim/DCSim: DCSim simulator release v0.3_ , https://doi.org/10.5281/ZENODO.8300961 (2023)