# A remark on MAKE – a Matrix Action Key Exchange

Chris Monico
Texas Tech University
<EMAIL_ADDRESS>Ayan Mahalanobis
IISER Pune
<EMAIL_ADDRESS>
##  1 Introduction

Rahman and Shpilrain [1] proposed a new key-exchange protocol MAKE based on
external semidirect product of groups. The purpose of this paper is to show
that the key exchange protocol is insecure. We were able to break their
challenge problem in under a second.

##  2 Description of MAKE

Let $\mathcal{G}$ and $\mathcal{H}$ be two semigroups of $k\times k$ matrices
over $\mathbb{F}_{p}$. The semigroup $\mathcal{G}$ is defined additively and
$\mathcal{H}$ is defined multiplicatively. We define a semidirect
$\mathcal{G}\ltimes\mathcal{H}$ such that
$(G_{1},H_{1})\cdot(G_{2},H_{2})=(H_{2}G_{1}H_{2}+G_{2},H_{1}H_{2})$, where
$G_{i}\in\mathcal{G}$ and $H_{i}\in\mathcal{H}$ for $i=1,2$. For more on the
description of the semidirect product used to make MAKE see [1].

Like all key exchange protocols, the purpose of MAKE is for Alice and Bob to
set up a common key for secure communications over an insecure channel. In the
case of MAKE this was achieved in the following way:

  1. (i)

Alice and Bob decide over an insecure channel that they are going to use two
invertible matrices $M$ and $H$ over $\mathbb{F}_{p}$ for some suitable prime
$p$. Here $M\in\mathcal{G}$ and $H\in\mathcal{H}$.

  2. (ii)

Alice chooses an integer $m$ and Bob an integer $n$. These integers are
private information.

  3. (iii)

Alice computes $(M,H)^{m}=(A,H^{m})$. The operation is the product in the
semidirect product defined above. She sends $A$ to Bob but keeps $H^{m}$
secret.

  4. (iv)

Bob computes $(M,H)^{n}=(B,H^{n})$ and sends $B$ to Alice and keeps $H^{n}$
secret.

  5. (v)

Alice on receiving $B$ computes the first component of
$(B,Q)\cdot(A,H^{m})=(H^{m}BH^{m}+A,QH^{m})$. The common key is
$H^{m}BH^{m}+A$. Note that $Q$ is neither known nor needed in this case and
the key exchange is successful without an explicit knowledge of $Q$.

  6. (vi)

Bob on receiving $A$ computes the first component of
$(A,Q)\cdot(B,H^{n})=(H^{n}AH^{n}+B,QH^{n})$. The key is $H^{n}AH^{n}+B$. Note
as before $Q$ is not known.

To check that this key exchange is successful one has to check if

$H^{m}BH^{m}+A=H^{n}AH^{n}+B$

is true. This follows from $A=H^{m-1}MH^{m-1}+H^{m-2}MH^{m-2}+\ldots+HMH+M$
and $B=H^{n-1}MH^{n-1}+H^{n-2}MH^{n-2}+\ldots+HMH+M$, and the common key is

$K=H^{m+n-1}MH^{m+n-1}+H^{m+n-2}MH^{m+n-2}+\ldots+HMH+M.$

##  3 An attack on MAKE

Our attack is based on the following lemma, illustrating that recovery of the
private parameters is not necessary for an eavesdropper to obtain the shared
secret key.

######  Lemma 1.

Let $M,H,A,B,m,n$ be as above. Suppose $R$ and $S$ are matrices which commute
with $H$ and satisfy

$RMS=HAH+M-A.$

Then Alice and Bob’s shared secret key $K$ satisfies $RBS+A=K$.

###### Proof.

We simply compute

$\displaystyle RBS+A$ $\displaystyle=$ $\displaystyle
R\left(\sum_{i=0}^{n-1}H^{i}MH^{i}\right)S+A$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{n-1}H^{i}RMSH^{i}+A$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{n-1}H^{i}(HAH+M-A)H^{i}+A$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}H^{i}AH^{i}+\sum_{i=0}^{n-1}H^{i}MH^{i}-\sum_{i=0}^{n-1}H^{i}AH^{i}+A$
$\displaystyle=$ $\displaystyle
H^{n}AH^{n}+\sum_{i=0}^{n-1}H^{i}MH^{i}=\sum_{i=0}^{m+n-1}H^{i}MH^{i}=K.$

∎

Our approach to finding such matrices $R,S$ is as follows. We will find
polynomials $f,g\in\mathbb{F}_{p}[t]$ for which

$f(H)M=\left(HAH+M-A\right)g(H),$ (1)

with $g(H)$ invertible. If the multiplicative order of $H$ is $\nu$, then
$f(t)=t^{m}$ and $g(t)=t^{\nu-m}$ satisfy (1), so such polynomials necessarily
exist. In that case, it follows that $R=f(H)$ and $S=g(H)^{-1}$ satisfy the
hypothesis of Lemma 1.

To find $f,g$ satisfying (1), first note that by the Cayley-Hamilton Theorem,
we may assume they each have degree at most $k-1$. Set $Z=HAH+M-A$, and is
then sufficient to solve

$\displaystyle
f_{0}M+f_{1}HM+\cdots+f_{k-1}H^{k-1}M=g_{0}Z+g_{1}ZH+\cdots+g_{k-1}ZH^{k-1}.$
(2)

This is a homogeneous linear system of $k^{2}$ equations in the $2k$ unknowns
$f_{0},\ldots,f_{k-1},g_{0},\ldots,g_{k-1}$. We use standard techniques to
produce a basis for the subspace of solutions in $\mathbb{F}_{p}^{2k}$. Then
simply choose nonzero solutions at random until one is found for which the
resulting $g(H)=g_{0}I+g_{1}H+\cdots+g_{k-1}H^{k-1}$ is invertible. We
implemented this attack based on the Python code provided by the authors of
[1] for the case of $k=3$. After performing 200 experiments for each of
several different primes ranging from $p=17$ to their 2000-bit prime number,
we did not encounter a single case where more than one random choice was
necessary. It’s not hard to show that if the dimension of the solution space
is 1, then every nonzero solution yields an invertible $g$; but we do not know
if that dimension is always 1.

##  4 Analysis of the attack

The attack presented here has just four steps:

  1. 1.

Compute $Z$.

  2. 2.

Find a basis for the solutions to (2),

  3. 3.

Choose a nonzero solution to (2) and determine $R=f(H)$ and $S=g(H)^{-1}$.

  4. 4.

Compute $K=RBS+A$.

Finding $Z$ requires 2 matrix multiplications and 3 additions, for a total of
$\mathcal{O}(k^{3})$ arithmetic operations in $\mathbb{F}_{p}$.

To compute a basis for the solutions to (2), one explicitly computes
$M,HM,H^{2}M,\ldots,H^{k-1}M$ using $k-1$ matrix multiplications, and then
$Z,ZH,ZH^{2},\ldots,ZH^{k-1}$ using another $k-1$ matrix multiplications, for
a total of $\mathcal{O}(k^{4})$ arithmetic operations in $\mathbb{F}_{p}$.
These $2k$ matrices are ‘flattened’ into row vectors to form the
$2k\times(k^{2}+2k)$ matrix consisting of the $2k\times k^{2}$ left submatrix
formed by those row vectors, augmented with the $2k\times 2k$ identity matrix.
This matrix is row-reduced, using
$\mathcal{O}(4k^{2}(k^{2}+2k))=\mathcal{O}(k^{4})$ arithmetic operations in
$\mathbb{F}_{p}$. A basis for the kernel is then found on the right-hand side
of the reduced matrix, as the row-vectors whose corresponding left-hand sides
are zero. Computing and inverting $g(H)$ requires $\mathcal{O}(k^{3})$ more
arithmetic operations in $\mathbb{F}_{p}$. Supposing that this must be
repeated $t$ times before obtaining an invertible $g$, we have a total of
$\mathcal{O}(tk^{4})$ arithmetic operations in $\mathbb{F}_{p}$, which
dominates the time required to compute $K$, so this is the total number of
$\mathbb{F}_{p}$ arithmetic operations required, with the most expensive
operation being inversion in $\mathbb{F}_{p}$. Therefore the process requires
no more than $\mathcal{O}(tk^{4}log^{3}{p})$ bit operations.

We performed 200 experiments with $k=3$ for each of the four primes
$p=17,19,135257,p_{2000}$, where $p_{2000}$ is the 2000-bit prime used by the
authors of [1] in their sample code. We did not encounter a case in which
$t>1$. So we experimentally conclude that the entire attack requires
$\mathcal{O}(k^{4}\log^{3}{p})$ bit operations, which is less than the cube of
the input size, and hence polynomial-time. On a single core of an i7-5557U
3.10 GHz processor, the time required to solve the largest $p_{2000}$
instances was about 0.03 seconds.

## References

  * [1] Nael Rahman and Vladimir Shpilrain.  MAKE: a Matrix Action Key Exchange.  Technical report, arXiv:2009.00716, 2020.


Simple analytic model for radiotherapeutic X-ray induced acoustic signal as a
function of absorption parameters

Canberk Sanli1,*, Esra Aytac Kipergil1, Mehmet Burcin Unlu1,2,

1 Department of Physics, Bogazici University, 34342 Bebek, Istanbul, Turkey
2 Center for Life Sciences and Technologies, Bogazici University, 34342 Bebek,
Istanbul, Turkey

*<EMAIL_ADDRESS>
## Abstract

Previous studies showed that it is possible make image reconstruction based on
the dose dependence of the therapeutic XA (X-ray induced acoustic signal)
amplitude which is then used to make dose mapping. We aimed to bring further
explicit parametrization for the acoustic signal in terms of the absorption
parameters, since this would mean encoding more information regarding the
absorption process to XA signals. The first step is to obtain pressure
waveform due to a point dose absorption by solving the thermo-acoustic
equation governing the heat absorption-pressure induction process based on the
analytic integration technique. Then, clinically relevant XA signal profile at
the detection point is obtained by generalizing point-dose-gradient induced
acoustic signal to surface-dose-gradient of a uniform spherical 3D dose
distribution based on the reciprocity principle for pressure waves in fluid
media. Therapeutic XA signal induced from the surface of the uniform spherical
dose distribution due to X-ray irradiation onto $5x5$ $cm^{2}$ field of the
water surface by $1\leavevmode\nobreak\ \mu s$ pulses delivering $1.7$
mGy/pulse is simulated in time and frequency domain. XA waves obtained in
previous empirical studies are simulated and compared by means of shape and
relative amplitude. Considering the previous studies on this subject, we
believe that the significance of this study is the foundation of a novel and
self-contained analytic approach to simulate the therapeutic X-ray acoustic
waves based on the physical parametrization of the energy transfer process.
This not only provides a better understanding of the physical phenomena
underlying the medical technique in terms of the medically relevant parameters
such as field size, pulse duration, absorbed dose per pulse etc. together with
the physical assumptions used to obtain a solution to the photo-acoustic
equation, but also brings consistent simulation results with previous
experimental and k-Wave results.

##  I Introduction

Pulsed X-ray beams produced in medical linacs (MeV range) or in X-ray tubes
(keV range) induce ultrasonic signals Bowen et al. (1991); Xiang et al. (2016)
from which tomographic reconstruction of the irradiated region, called as
X-ray induced Computed Tomography (XACT), is possible Xiang et al. (2013,
2016); Hickling et al. (2014a). Hence, linac generated acoustic signals can be
used for in vivo dose monitoring during radiotherapy provided the higher
spatial resolutions are achieved Xiang et al. (2013), and the latter can be
used to achieve deeper imaging than the current laser-based photoacoustic
imaging techniques supply provided that low dose absorption does not pose a
clinical risk Xiang et al. (2016). For both purposes, one needs to know the
dependence of the acoustic wave profile on the absorption parameters. The
linear dependence of pressure amplitude on the absorbed energy is already
known from the general photoacoustic theory Wang and Wu (2012) and
qualitatively verified through XACT images Hickling et al. (2016) by comparing
XACT and ion chamber images. The dependence on the pulse duration is also
qualitatively indicated by Xiang et al. (2016).

Besides the plausibility of implementing XACT in clinical practice as X-ray
imaging, X-ray induced acoustic (XA) waves can also be used as a real-time in
vivo dosimetry tool in radiotherapy. Since the X-ray sound solely depends on
the X-ray absorption, any microscopic change in tissue properties due to
radiotherapeutic beams can in principle be determined from the induced
acoustic signals Xiang et al. (2016) if a complete description of the X-ray
absorption regarding X-ray beam and tissue parameters is achieved. This would
be an invaluable tool as it can increase the efficiency of radiotherapy
treatments by monitoring the instantaneous treatment respond which may
necessitate unplanned interventions. The feasibility of using XA waves as a
dosimetry tool is emphasized in previous medical physics studies Xiang et al.
(2013); Hickling et al. (2014a); Xiang et al. (2016), and presciently by
Askariyan Askaryan (1957) when he first mentioned the generation of acoustic
signals due to the passage of ionizing radiation through matter.

Hickling et al. (2016) investigated this possibility by assuming a spatial
profile for the heating function $H(\mathbf{r})$ which describes the
instantaneous localized heating due to radiation exposure and defined through
a ‘heat defect function’ $k(\mathbf{r})$ in which all external absorption
parameters such as pulse duration of the X-ray beam and absorption cross-
section are implicit. Based on this formalism, they composed XACT images of
which intensity profile formed by the relative pressure amplitudes of the
ultrasound waves is used to obtain the dosimetric information which is then
verified by separate measurements performed in an ion chamber. Consequently,
the dose distribution is simulated with an accuracy depending on that in
determining the heat defect, physical density, Grüneisen coefficient, and
properties of detection apparatus.

This approach of using XACT images formed by the relative pressure amplitudes
which are determined by the heating function defined above has two main
theoretical deficiencies if its implementation for real-time in vivo imaging
in radiotherapy is considered. Firstly, since XACT images are based only on
the relative pressure amplitudes, all the remaining information about the
pressure signal is completely lost. This includes width and exact shape of the
pressure profile. Secondly, the above definition for the heating function
yields an absorption profile which does not give an explicit connection
between X-ray beams dependent absorption parameters such as pulse duration and
absorption cross section. Hence, one cannot relate the above-mentioned XACT
images with these parameters which are important in dose absorption during
radiotherapy Xiang et al. (2016).

These deficiencies do not cause any harm if XA waves are only to be used for
imaging structural heterogeneities in matter such as lead blocks of different
shapes in water Hickling et al. (2014b, 2016), in chicken breast Xiang et al.
(2013), or tiny piece of chicken bone in water phantom Xiang et al. (2016).
However, there exists no such contrast in shape and structure between normal
and tumor tissue. Indeed, given that they are both treated as a water medium
for most purposes, there exists no a priori reason why the amplitudes of the
XA signals should significantly differ between tumor and normal tissue. Hence,
XACT, as it is proposed, cannot be used to determine whether a detected XA
signal is induced by the dose on a normal or tumor tissue which would be the
main reason why one might perform a dose mapping during or after radiotherapy.

The aim of measuring induced acoustic signal is obtaining information
concerning the absorption process, and the complete description of the X-ray
absorption by deciphering the acoustic signal is possible, in principle. This
means provided that accurate theoretical model with sufficient parameters is
proposed, one can distinguish between XA signals due to different amount of
dose absorption and different characteristics of medium where absorption takes
place. Hence, an effort should be made to optimally model the dependence of
pressure profile on the absorption process regarding the tissue and X-ray beam
properties. Previous studies showed that it is possible make image
reconstruction based on the dose dependence of pressure amplitude which is
then used to make dose mapping. We aimed to bring further (than dose/amplitude
dependence) explicit parametrization for acoustic signal in terms of the
absorption parameters, since this would mean encoding more information
regarding the absorption process to XA signals. For this purpose, we propose a
novel analytic model for induced acoustic signals from the surface of a
uniform spherical dose distribution due to radiotherapeutic X-ray exposure
inside the tissue, and observed at an arbitrary position outside the
distribution inside the same medium. The first step is to obtain pressure
waveform due to a point dose absorption by solving the thermo-acoustic
equation governing the heat absorption-pressure induction process based on the
analytic integration technique. Then, clinically relevant XA signal profile at
the detection point is obtained by generalizing point-dose-gradient induced
acoustic signal to surface-dose-gradient of a uniform spherical 3D dose
distribution based on the reciprocity principle for pressure waves in fluid
media Landau and Sykes . As a result, just like pressure amplitude mainly
changes with the absorbed dose, we showed that pulse duration and absorption
cross section also governs the waveform of the acoustic signal. It is an
interesting curiosity whether these explicit dependencies of XA signal can be
translated into the image reconstruction technique.

Thermo-acoustic model is first proposed to explain the induced acoustic
signals due to energy transfer during the passage of high energy particles
through matter; independently by Askariyan et al. (1979) and Bowen and Learned
(1979), and verified by Sulak et al. (1979). The solution to the thermo-
acoustic equation is yielded by considering the deposited heat power per unit
volume as a fictitious charge distribution and the observed pressure as
‘retarded potential’ Askariyan et al. (1979); Sulak et al. (1979); Griffiths
(2005). Then, by assuming a Gaussian distribution for instantaneous heat
deposition, dependencies between pressure profile and particle cascade
parameters are obtained Learned (1979) by direct integration. Only physical
assumptions used to simplify the integration giving the point-source pressure
wave solution are instantaneous heat deposition $\dot{T}˙=T(r)\delta(t)$ due
to X-ray exposure where $T(r,t)$ is the temperature rise at position $\vec{r}$
on the treatment plane, and sufficiently wide-band transducers so that
absorption is frequency independent Askariyan et al. (1979). Validity of this
technique in medical linac acoustic cases is already justified in medical
physics literature Baily (1992). Former physical assumption being true for any
ionizing radiation can be elucidated by showing that the thermal relaxation
time determined by the minimum dimensions of the area of heat evolution is
essentially larger than the characteristic time of action of penetrating
radiation determined by the pulse duration Lyamshev (2004), so that heat
conduction can be neglected Wang et al. (2003). As shown in the following
section, this relation due to the implicit assumption of no sheer wave
generation in the thermo-acoustic equation also serves as a bridge in
expressing the absorption cross section in terms of a width parameter of a
Gaussian curve which is a fit to the $5x5\leavevmode\nobreak\ cm^{2}$ square
field on the treatment surface and forms the area of heat evolution. On the
other hand, detection by immersion ultrasound transducer with a sufficiently
large bandwidth so that frequency dependence of the absorption can be ignored
is an priori assumption of our model. This is a reliable assumption
considering the sharp peak of the therapeutic XA signal at a central frequency
Hickling et al. (2014b).

Although there have been many analytic simulations of the medical proton-
acoustic emission after the confirmation of the thermo-acoustic model for
treatment beams Baily (1992), systematic study of dependency between physical
absorption parameters and acoustic wave profile parameters is lacking Jones et
al. (2014). Analogously, Hickling et al. (2014b) introduced a simulation
workflow combining dose simulations and the induced acoustic wave transport
simulation to form XACT images from which dosimetric information is extracted.
This approach, however, does not reveal the explicit dependence between
pressure profile and other important absorption parameters pulse duration and
absorption cross section Xiang et al. (2016).

The method of solving thermo-acoustic equation is to assume a Gaussian
distribution for the instantaneously deposited heat profile of which volume
integral is then evaluated over the region of heat deposition at the retarded
time, analogously to the derivation of analytic dependence of acoustic waves
on absorption parameters from the thermo-acoustic model of high energy
particles Learned (1979); Askariyan et al. (1979). Analytic expression for the
acoustic pressure profile generated by a point dose deposition inside the
sample is then obtained as a function of the absorbed dose at that point, and
absorption parameters pulse duration and absorption cross section. Dose
gradient rate inducing acoustic pulses are determined by both beam energy and
pulse repetition frequency, whereas pulse duration is solely determined by the
linac operation Karzmark and Pering (1973). On the other hand, temperature
rise at any point of the irradiated medium depends on the local absorption
probability governed by the absorption cross section $\sigma_{abs}$ which is
given by the product of mass energy attenuation coefficient with the atomic
number density of the material. We approximate this temperature rise profile
by a Gaussian curve of width parameter $\sigma$ which is fit to the square
field size on the treatment surface and determines the minimum dimension of
the area of heat evolution which is equal to the inverse of the mass energy
attenuation coefficient for X-ray absorption Lyamshev (2004). Hence, given a
uniform medium, absorption cross section for X-rays only depends on the field
size on the surface which on the other hand defines the width parameter of
Gaussianly distributed heat absorption profile in our model. In other words,
given a uniform medium, absorption cross section only depends on the width
parameter of Gaussianly distributed heat absorption profile in our model.

In a clinically applicable case the acoustic signal reaching the detector is a
superposition of this point-wise generated waves in accordance with the 3D
dose distribution due to X-ray exposure, and the detector is located at an
arbitrary position outside this distribution. Hence, we generalized the
analytic expression for point-dose-gradient induced acoustic signal during
X-ray radiotherapy to surface-dose-gradient of a uniform spherical 3D dose
distribution based on the reciprocity principle, and ignoring the attenuation
of the acoustic signal due to small attenuation coefficient, exponential decay
and close-range detection. This generalization enables the consistent
comparison of the simulated acoustic waveforms based on our analytic model
with the plausible clinical scenario where isodose curve is approximately
spherical due to depth-dose curve for X-rays Townson and Zavgorodni (2014) and
induces acoustic signals as previously observed and simulated by Hickling et
al. (2016). Furthermore, dependence of the XA signal on the dose rate,
absorption cross section, and pulse duration are quantitatively investigated
by simulations in the therapeutic range of these parameters.

The approach taken here can be used for any radiation induced acoustic wave
mechanism in general, given the expression for the initial pressure rise
depending on the local energy absorption properties, particularly for proton
induced acoustic waves of which the initial pressure rise is similarly given
Ahmad et al. (2015).

## Materials and Methods

For a point source of heat which is located at $\vec{r}$ and described by a
macroscopic heat deposition function $H(\vec{r},t)$ associated with the
temperature rise $T(\vec{r},t)$ at any point $\vec{r}$ on the calculation
plane which is caused by the dose gradient around a point $\eta_{0}$ taken as
the origin of our coordinate system, the generation and propagation of XA wave
are given by the photoacoustic equation Wang and Wu (2012); Wang (2008); Morse
and Ingard (1968)

$\displaystyle(\nabla^{2}-\frac{1}{v_{s}^{2}}\frac{\partial^{2}}{\partial
t^{2}})p(\vec{r},t)$ $\displaystyle=$ $\displaystyle-\frac{\beta}{\kappa
v_{s}^{2}}\frac{\partial^{2}T}{\partial t^{2}}$ (1)
$\displaystyle(\nabla^{2}-\frac{1}{v_{s}^{2}}\frac{\partial^{2}}{\partial
t^{2}})p(\vec{r},t)$ $\displaystyle=$ $\displaystyle-\frac{\beta\rho
C_{v}}{C_{p}}\frac{\partial^{2}T}{\partial t^{2}}$ (2)

of which solution $p(\vec{r},t)$, i.e. the pressure at any point $\vec{r}$ and
at time $t$, is most simply expressed in terms of the volume integral over the
region of absorption (heat deposition), i.e. space $V$ which is formed by a
point source of heat located at $\vec{r^{\prime}}$ and bounded by the stress
confinement sphere $S_{r}^{R}=v_{s}\tau$ where $v_{s}$ is the speed of sound
in the medium and $\tau$ is the pulse duration, at the retarded time
$t-{|\vec{r}}-{\vec{r^{\prime}}|}/v_{s}$ Sulak et al. (1979); Baily (1992);
Terunuma et al. (2007); Kruger et al. (1995)

$p(\vec{r},t)=\frac{\beta\rho C_{v}}{4\pi
C_{p}}\int_{V}\frac{dV^{\prime}}{{|{\vec{r}}-{\vec{r^{\prime}}}}|}\frac{\partial^{2}}{\partial
t^{2}}T{\bigg{(}}{\vec{r^{\prime}}},t-\frac{{|\vec{r}}-{\vec{r^{\prime}}|}}{v_{s}}{\bigg{)}}$
(3)

which can be further reduced to the integration over the spherical surface of
radius $R=v_{s}\tau$ with its center at $\vec{r}$ Askariyan et al. (1979)

$p(\vec{r},t)=\frac{\beta\rho C_{v}}{4\pi
C_{p}}v_{s}^{2}\frac{\partial}{\partial
R}\int_{S_{r}^{R}}dS\frac{T(\vec{r},t)}{R}$ (4)

by assuming instantaneous heat deposition $\dot{T}˙=T(r)\delta(t)$ due to
X-ray exposure, and sufficiently wide-band transducers so that absorption is
frequency independent, where $\beta$ is thermal coefficient of volume
expansion, $\rho$ is the mass density of soft tissue, and $C_{p}$ and $C_{v}$
denote the specific heat capacities at constant pressure and volume which
together define the isothermal compressibility $\kappa=C_{p}/(\rho
v_{s}^{2}C_{v})$.

Figure 1: Treatment Field and Integration Region

Here, the temperature rise $T(\vec{r},t)$ is related to the heating function
$H(\vec{r},t)$ defined as the thermal energy converted per unit volume per
unit time by Wang and Wu (2012)

$\rho C_{v}\frac{\partial T(\vec{r},t)}{\partial t}=H(\vec{r},t)$ (5)

which is expressed in terms of a product of spatial and temporal Gaussian
profiles as

$H(\vec{r},t)=\frac{p_{0}C_{p}}{\beta
v_{s}^{2}}e^{-r^{2}/2\sigma^{2}}\frac{1}{\sqrt{2\pi\tau^{2}}}e^{-t^{2}/2\tau^{2}},$
(6)

on the macroscopic treatment surface taken as the $5x5\leavevmode\nobreak\
cm^{2}$ square field, in order to reveal the dependence of pressure profile on
the X-ray absorption parameters: pulse duration-$\tau$, and absorption cross
section-$\sigma_{abs}$ which is defined by the width parameter $\sigma$ of the
Gaussian fit to the treatment surface.

Integrating Eq (6) gives

$\rho C_{v}T(\vec{r},t)=\frac{p_{0}C_{p}}{\beta
v_{s}^{2}}\frac{e^{-r^{2}/2\sigma^{2}}}{\sqrt{2\pi\tau^{2}}}\int_{0}^{t}e^{-t^{\prime
2}/2\tau^{2}}dt^{\prime}$ (7)

where

$\int_{0}^{t}e^{-t^{\prime
2}/2\tau^{2}}dt^{\prime}=\tau\sqrt{\frac{\pi}{2}}\mathbf{erf}(t/\sqrt{2}\tau)$
(8)

and hence

$\rho C_{v}T(\vec{r},t)=\frac{p_{0}C_{p}}{\beta
v_{s}^{2}}\frac{e^{-r^{2}/2\sigma^{2}}}{2}\mathbf{erf}(t/\sqrt{2}\tau)$ (9)

Here, we note that $T(\vec{r},t)$ is not exactly constant on the integration
surface $S_{r}^{R}$ centered at r, with radius $R=v_{s}\tau$, as in Fig 1.
Since $R$ is defined by the region of stress confinement that is defined by
the pulse duration, and hence is small whereas $r$ is always of the
macroscopic size; $T(\vec{r},t)$ is assumed to be constant and moved outside
the integration in Eq (4) so that after substituting (9) into (4) we get

$\displaystyle p(\vec{r},t)$ $\displaystyle=$
$\displaystyle\frac{p_{0}e^{-r^{2}/2\sigma^{2}}}{8\pi}\mathbf{erf}(t/\sqrt{2}\tau)\frac{\partial}{\partial
R}\int_{S_{r}^{R}}R^{2}(1/R)sin\theta d\theta d\phi$ (10) $\displaystyle
p(\vec{r},t)$ $\displaystyle=$
$\displaystyle\frac{p_{0}}{2}e^{-r^{2}/2\sigma^{2}}\mathbf{erf}(t/\sqrt{2}\tau)$
(11)

To sum up, the idea is to consider the macroscopic temperature rise
distribution as a Gaussian curve lasting for a small time of a pulse duration
$\tau$ defining at any point $\vec{r}$ a small region of stress confinement
$S_{r}^{R}$ of radius $R=v_{s}\tau$ over whose surface $T(r,t)$ remains
constant. Thus, Eq (11) which can also be seen as the special case of the
previously established photoacoustic deconvolution algorithm Wang et al.
(2004) with Gaussian deposition of heat Learned (1979) gives the pressure rise
at point $r$ due to initial pressure rise $p_{0}$ at the origin. This defines
a pressure profile (field) $p(r,t)$ due to a single point heat source over the
macroscopic treatment region restricted by the field size
$5x5\leavevmode\nobreak\ cm^{2}$. We neglect the attenuation of the induced
ultrasound due to small attenuation coefficient, exponential decay and close-
range detection to conclude that a transducer located nearby this region would
observe the same profile $p(r,t)$. We will generalize this pressure profile
induced by a point-dose-gradient to the waveform $p^{3D}(\vec{\eta},t)$
induced from a dose gradient on the surface of a uniform dose distribution as
a function of the distance to the surface of the sphere from the transducer
location $\vec{\eta}$.

For an X-ray induced acoustic wave (XA), the initial pressure rise is given by
Hickling et al. (2014b)

$p_{0}(\vec{r})=\Gamma\rho\tilde{D}(\vec{r})$ (12)

where $\tilde{D}$ is the absorbed dose per pulse at the position $\vec{r}$ on
the treatment surface, and $\Gamma$ being the Gruneisen parameter describing
conversion efficiency between deposited heat energy and pressure Hickling et
al. (2016) is defined in terms of thermal coefficient of volume expansion
$\alpha$ and specific heat capacity at constant pressure $C_{p}$ as
$\Gamma=\frac{\beta}{C_{p}}v_{s}^{2}$ where $v_{s}$ is the speed of sound in
the medium, and taken to be uniform; $\Gamma=0.15$, since detection is assumed
to be realized in a homogeneous water tank at room temperature. Hence,
acoustic pressure profile of XA wave at ($\vec{r}$,t) due to point heat source
located at $\vec{r}_{0}$ is given by

$\displaystyle p(\vec{r},t)$ $\displaystyle=$
$\displaystyle\frac{\Gamma\tilde{D}(\vec{r}_{0})\rho}{2}\leavevmode\nobreak\
e^{-r^{2}/2\sigma^{2}}\leavevmode\nobreak\ \mathbf{erf}(t/\sqrt{2}\tau)$ (13)

where $\tilde{D}(\vec{r}_{0})=\tilde{D}$ is the absorbed dose per pulse at
point $\vec{r}_{0}$ on the treatment surface and is the first order
approximation to the nearby dose in the region of interest, i.e.
$D(\vec{r})\simeq D(\vec{r}_{0})$ for ${r}<5\leavevmode\nobreak\ cm$, and
determined by the medical range of dose rates and effective pulse repetition
frequencies at which linacs are operated. To avoid any confusion concerning
the range of the parameters, we again emphasize that $t$ defines a microscopic
time related with the stress confinement whereas $r$ defines the macroscopic
region characterized by a Gaussian of width $\sigma$ given by the field-size
on the surface, and they together form the temperature profile in space and
time of which surface integral at each point over the region of stress
confinement as in Eq (4) forms the pressure profile on the treatment surface
in the course of several pulse durations.

In a realistic case however the acoustic signal reaching the detector would be
the superposition of our point-wise generated waves given by Eq (13) in
accordance with the 3D dose distribution due to X-ray exposure, and the
detector would be in an arbitrary position from the dose gradient that induces
XA signal. To achieve this, we first assumed 3D spherical uniform dose
distribution inside the sample in a region of $l_{min}^{3}\leavevmode\nobreak\
cm^{3}$ where $l_{min}$ is the minimum dimension of the area of heat evolution
being consistent with the relative dose-depth curve for X-rays Townson and
Zavgorodni (2014); Podgorsak (2005) together with the assumption of isotropic
media, and determined by the radius of the radiation spot $a/2$ for X-ray
irradiation by $l_{min}\approx a/2\approx\sigma$ Lyamshev (2004) where $a$ is
the field size, as will be explained. Pressure waves are then induced due to
dose gradient on the surface $S$ of this spherical region. Our aim is to
calculate the pressure measured by a transducer located nearby the sphere. Let
us put the transducer first at the center of the sphere. According to the
reciprocity principle, pressure at this center due to all point sources (dose
per pulse) located on the sphere is equal to the total pressure on the
spherical surface due to point source (dose per pulse) located at the center.
This is equivalent to assuming that point sources on the sphere is distributed
homogeneously. Hence, to find the pressure at the center of the spherical dose
distribution, $p^{3D}(\vec{r_{0}},t)$, we equivalently find the average
pressure field due to a point source Eq. (12) on the surface $S$.

$\displaystyle p^{3D}(\vec{r_{0}},t)$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi(a/2)^{2}}\int_{S}\frac{\Gamma\tilde{D}(\vec{r_{0}})\rho}{2}\leavevmode\nobreak\
e^{-\eta^{2}/2\sigma^{2}}\leavevmode\nobreak\ \mathbf{erf}(t/\sqrt{2}\tau)da$
(14) $\displaystyle p^{3D}(\vec{r_{0}},t)$ $\displaystyle=$
$\displaystyle\frac{\Gamma\tilde{D}\rho}{2}\leavevmode\nobreak\
\frac{\eta^{2}}{\sigma^{2}}\leavevmode\nobreak\
e^{-\eta^{2}/2\sigma^{2}}\leavevmode\nobreak\ \mathbf{erf}(t/\sqrt{2}\tau)$
(15)

where $\vec{\eta}$ is the distance to the surface of the sphere from the
transducer location, and equal to the radius of the sphere for the case where
the transducer is located at the center.

Second approximation is that ultrasound attenuation is negligible due to small
attenuation coefficient, exponential decay and close-range detection. Hence,
as we move the transducer radially outside the sphere, at a radial distance
$\vec{\eta}$ from the surface, it will measure the signal given by

$p^{3D}(\vec{\eta},t)=\frac{\Gamma\tilde{D}\rho}{2}\leavevmode\nobreak\
\frac{\eta^{2}}{\sigma^{2}}\leavevmode\nobreak\
e^{-\eta^{2}/2\sigma^{2}}\leavevmode\nobreak\ \mathbf{erf}(t/\sqrt{2}\tau)$
(16)

Eq (14) is a general expression for the pressure profile induced by the dose
gradient at the surface $S$ of a 3D uniform dose distribution with a center at
$r_{0}$ taken as the detection point. For a uniform spherical dose
distribution, it yields Eq (16) as the pressure profile measured by the
transducer at a distance $\vec{\eta}$ from the surface of the sphere. We note
that $\eta$ takes values between $d$ and $2l_{min}+d$ where $d$ is the closest
distance to the sphere from the transducer and $2l_{min}=a$ is the field size
and forms pressure profile at the transducer location as given below. This can
be consistently compared with the previous X-ray induced acoustic studies,
particularly with the study by Hickling et al. (2016) where both experimental
and k-wave simulated acoustic signals reaching from boundaries of various
uniform 3D dose distributions in pure water due to different collimation of
the X-ray field.

Dose gradient per pulse $\tilde{D}$ in Eq (16) which is the main parameter of
interest in X-ray acoustics can be determined via two approaches. One
alternative is to use clinical dose values per fraction which are required in
X-ray radiotherapy to be in between $1.7$ and $2.5$ Gy Maciszewski and Scharf
(2004) and are verified through various experimental techniques Podgorsak
(2005); Mayles et al. (2007); Constantinou et al. (1982); Quach et al. (2000);
Zhu et al. (1995); Butson et al. (1996); O’Keeffe et al. (2016), and also by
Monte-Carlo simulations and analytical modeling Rogers et al. (1995); Ding
(2002); Jagetic and Newhauser (2015); Taddei et al. (2013). This determines a
dose gradient at the surface of the dose distribution due to the difference
between the amount of dose in the treatment and the nearby regions. This
difference is approximately 2 Gy fraction Jagetic and Newhauser (2015) and
corresponds to $\tilde{D}=1.7$ mGy/pulse by taking the effective pulse
repetition frequency as 60 Hz Karzmark and Pering (1973); O’Keeffe et al.
(2016) and macro-pulse duration to deliver 2 Gy fraction as 20 seconds
O’Keeffe et al. (2016). Alternatively, one can use the pulse repetition
frequency value directly with the dose rate taken as 6 Gy/min in the
operational range at which medical linacs are operated in radiotherapy
Maciszewski and Scharf (2004) to get the same dose-per-pulse value. Either way
we deduce that radiotherapeutic dose gradient per pulse on the sphere can vary
in the the interval 1 to 2 mGy/pulse as also taken in the previous therapeutic
XA studies Xiang et al. (2013); Hickling et al. (2014b, 2016) and we simulate
the acoustic signals accordingly.

Ultrasound-absorption dependence is investigated in terms of three parameters:
absorbed dose per pulse, pulse duration and the absorption cross section.
Range of pulse duration is taken in between 1-5 $\mu s$ as this is the
standard range in which medical linacs are operated since their first
appearance Karzmark and Pering (1973); O’Keeffe et al. (2016). Heat deposition
function (and hence the temperature rise) function $H(\vec{r},t)$ is
approximated as a Gaussian curve defined on all the calculation plane in Fig
1, as given by Eq 6. The width parameter $\sigma$ of this Gaussian is
calculated as

$\sigma=\frac{a}{2\sqrt{2ln2}}$ (17)

by fitting a Gaussian curve centered at the center of the square and of which
full width at half maximum intersects the edge of the square field of size
’a’. On the other hand, since absorption cross section is given by
$\sigma_{abs}=\mu A/(\rho N_{A})$, where $\rho$ is the mass density, $\mu$ is
the X-ray absorption coefficient, and $N_{A}$ and $A$ are Avogadro number and
atomic number, respectively, it is parametrized only by attenuation
coefficient $\mu$, or equivalently by the minimum dimension of the area of
heat evaluation $l_{min}$ given by $l_{min}=\mu^{-1}\approx a/2$ for X-rays
Lyamshev (2004)for a given medium, where $a/2$ is the radius of the radiation
spot at the surface depending on the size of the treatment field which is
taken as $5x5\leavevmode\nobreak\ cm^{2}$ in this study. The existence of this
minimum dimension which is at the order of the field size (or $\sigma$) on the
surface is due to the implicit assumption in the photo-acoustic equation that
pulse duration is small enough to neglect the viscous effects so that we have
$l_{min}>>1/(\tau\nu)^{1/2}$, i.e. the minimum dimensions of the area of heat
evolution (thermal source) must be essentially larger than the depth of wave
penetration which is characterized by $1/(\tau\nu)^{1/2}$ where $\tau$ is
pulse duration and $\nu$ is kinematic coefficient of viscosity. Therefore, in
X-ray radiotherapy, assuming a uniform medium, absorption cross section is
defined solely by the field size, and hence by the Gaussian width parameter
$\sigma$ of the heating function. This means that variation in the $\sigma$
can be interpreted in two ways. That could be taken as a change in the field
size, or a change in the medium, once one of them is fixed. In both cases, the
effect can be observed in the pressure waveform and simulated by our simple
analytic model. Corresponding pressure profiles as a function of dose gradient
on the surface of a uniform spherical dose distribution per pulse, pulse
duration, and absorption cross section as given by Eq (16) are plotted using
MATLAB, Mathworks Inc. Based on these profiles, the central frequency $f_{0}$
of the signal is determined by applying Fast Fourier Transform (FFT) from
which the noise level of a generic transducer with uniform efficiency
$\eta(f)\approx\eta(f_{0})=0.5$ and with size $A=3\leavevmode\nobreak\ cm^{2}$
to detect this signal is calculated by using Winkler et al. (2013); Xiang et
al. (2016)

$NEP=\sqrt{k_{B}T\Bigg{[}1+\frac{F_{n}}{\eta(f_{0})}\Bigg{]}Z_{a}/A}\leavevmode\nobreak\
.\leavevmode\nobreak\ \sqrt{f_{0}}$ (18)

where $k_{B}$ is Boltzmann constant, $T$ is the absolute temperature,
$F_{n}=2$ is the noise factor of the amplifier, and $Z_{a}$ is the
characteristic acoustic impedance of the medium taken as water.

Similarly, one can use Eq (13) in order to plot acoustic signals induced by
proton beams used in proton radiotherapy as well, by directly implementing the
dose values at the Bragg peak Koehler et al. (1975); Larsson (1961); Clasie et
al. (2010), which is also analytically approximated by Bortfeld (1997).

## Results

Results contain two parts. First, we investigated the pressure signals induced
by different amounts of dose gradients per pulse $\tilde{D}$ on the surface of
a uniform spherical dose distribution. For this purpose, we fixed pulse
duration $\tau=1\leavevmode\nobreak\ \mu s$ and the width parameter
$\sigma=2.12$ cm, and by plotting the Eq (16) for different values of dose per
pulse values in the clinical range of 1 to 2 mGy/pulse with
$\eta=-5:(15/2)10^{-4}:10$ and $t=-10^{-5}:10^{-9}:10^{-5}$, we obtained XA
signal measured by a transducer at a distance $\eta$ from the surface of the
sphere. Secondly, for a fixed amount of clinical dose gradient per pulse of
1.7 mGy/pulse on the surface, the dependencies of the pressure signal on
$\tau$ and $\sigma$ are obtained by plotting Eq (16) for pulse duration values
of $1-5\leavevmode\nobreak\ \mu s$ at a fixed $\sigma=2.12$ cm and for
different width parameters 2.12, 1, 0.5 cm at a fixed
$\tau=1\leavevmode\nobreak\ \mu s$.

A therapeutic XA signal as being defined by the above mentioned parameters is
found to have two crests and troughs corresponding to the compression and
rarefaction regions separated by negative and positive dose gradients at the
two boundaries of the uniform dose distribution. Pressure waveform as being
observed by a transducer at a distance $\eta_{ct}=10$ cm from the center of
the sphere is simulated by Eq (16) in time and frequency domain as it is shown
in the Fig 2. Based on the power spectrum of the signal which is yielded by
applying FFT to the waveform in time domain and determines the central
frequency as $45.8$ kHz, $NEP$ is calculated from Eq (18) as $2.2$ mPa.
Further, acoustic signals induced by X-ray pulses delivering 1.7 mGy dose per
pulse with pulse durations of $1-3\leavevmode\nobreak\ \mu s$ at a fixed
absorption cross section of $2.12\leavevmode\nobreak\ cm$ and those with a
fixed pulse duration of $1\leavevmode\nobreak\ \mu s$ with different
absorption cross sections defined by different width parameters of 1.5, 1.7,
and 2.12 cm are obtained as given in Fig 3.

Figure 2: Therapeutic XA signal

Figure 3: Dependence of XA signal on Absorption Parameters

Qualification of our model is tested by simulating the XA waves obtained in
empirical studies of Xiang et al. (2013), and Hickling et al. (2014b, 2016).
Necessary parameters to simulate the acoustic signals based on Eq (16) are
given in the Table 1 where $D_{pp}$, $\tau$, and F.S. denote dose per pulse,
pulse duration, and treatment field size. Xiang et al. (2013) operated 10 MV
linac to induce X-ray pulses from a lead rod embedded in a chicken breast
tissue at a dose rate $30$ Gy/min and the pulse repetition frequency of $128$
pps from which the corresponding dose per pulse value is calculated as $3.9$
mGy. Hickling et al. (2014b) operated 6 MV linac at a pulse repetition of
$180$ Hz to induce acoustic waves from a lead block suspended in a water tank
in “clinically relevant radiotherapy dosimetry situations” wherefore
Maciszewski and Scharf (2004) dose per pulse value can be taken as being in
the range $1$ to $2$ mGy/pulse, and recently Hickling et al. (2016) observed
and simulated the XA signals induced from various dose distributions formed by
primary collimator (or multi-leaf collimator for non-standard fields) inside
pure water due to 4 $\mu s$ X-ray pulses delivering 1.11 mGy/pulse to make
image reconstruction. We calculated the corresponding width parameters
$\sigma$ of the Gaussian absorption profiles from the field size dimension by
using Eq (17). Consequently, we have simulated the corresponding induced
acoustic signals using our analytic model as shown in Fig 4. Validation of our
model is mainly realized by comparing XA signal in Fig 4(a) with the
corresponding empirical signal Hickling et al. (2016) by means of shape and
relative amplitude. Other experimental XA signals Xiang et al. (2013);
Hickling et al. (2014b) as being induced from the irradiation of the metal
blocks are not suitable for a consistent comparison. However, assuming that
the medium of irradiation will not significantly change the generic shape of
the therapeutic XA signal since it is restrained by the medical range of the
operational parameters in which those experiments are performed, qualitative
comparison of the waveforms is made. Among these studies, only Hickling et al.
(2014b) shows experimental therapeutic XA signal induced from the therapeutic
X-ray irradiation of a lead block in physical units of 0.2 Pa maximum value,
whereas in the other studies experimental pressure amplitudes are given in
arbitrary units. This value serves as a definite upper limit for the amplitude
of the XA signal induced from pure water that is simulated in this study since
density of lead is ten times larger than water.

Table 1: Parameters to Simulate Empirical Studies

| Xiang et al. (2013) | Hickling et al. (2014b) | Hickling et al. (2016)
---|---|---|---
Energy (MeV) | 10 | 6 | 10
$\mathbf{D_{pp}}$ (mGy/pulse) | 3.9 * | 1-2** | 1.11
$\mathbf{\tau}$ $\mathbf{(\mu s)}$ | $5$ | $4$ | $4$
F.S. $\mathbf{(cm^{2})}$ | 4x4 | 10x10 | 4x4
$\mathbf{\sigma}$ $\mathbf{(cm)}$ *** | 1.7 | 4.25 | 1.7

  * •

* calculated from Dose rate (Gy/min) and pulse repetition frequency.

  * •

** in clinical range

  * •

*** calculated from field size by using Eq (17).

Figure 4: Simulation of Previous XA Experiments. (a) Hickling et al. (2016)
with transducer located at different distances of $\eta_{ct}$ from the center
of the field (b) Xiang et al. (2013) with irradiation realized on a pure
water.

## Discussion and Conclusion

Therapeutic XA signal induced from the surface of the uniform spherical dose
distribution due to X-ray irradiation onto $5x5$ $cm^{2}$ field of the water
surface by $1\leavevmode\nobreak\ \mu s$ pulses delivering $1.7$ mGy/pulse is
simulated in time and frequency domain as given in Fig 2, and also with
different operation parameters in the medical range as given in Fig 3. The
linear dependence on dose rate and pressure amplitude is evident in Fig 2(a)
and in Fig 3(a), and the x-axis limits are deliberately fixed to reveal the
flattening effect of increasing pulse duration on XA signal in Fig 3(b) which
can be related to the decrease in the area of heat evolution due to increased
rate of thermal relaxation. Lastly, assuming that the same dose is delivered
by the pulses of the same length, the shape of the waveform undergoes just a
phase shift and a displacement when the absorption cross section which is
defined by the width parameter of the Gaussian fit to the treatment field is
changed as in Fig 3(c).

To be able to justify the above observations, we simulated the proposed
analytic model Eq (16) by imposing the parameters used in previous XA studies.
Sanity check of our model is realized through comparisons in waveforms,
amplitudes, and power spectrum. Simulated therapeutic XA signal from the dose
gradient at the treatment edges inside water is given in Fig 4(a) as observed
from different transducer locations, and it has the same waveform with the
corresponding empirical signal Hickling et al. (2016) in which rarefaction
preceding compression and compression preceding rarefaction regions are
identified by crest-through and through-crest portions of the waveform which
are caused by the positive and negative dose gradients at the two boundaries.
Hence, we see that relative amplitude between a crest and trough gives the
dose gradient information at a boundary. Experimental therapeutic XA waveforms
induced from lead blocks Xiang et al. (2013); Hickling et al. (2014b) are also
found to be the similar in shape to the analytically simulated signal.
Secondly, we note that the amplitude of analytically simulated therapeutic XA
signal induced from pure water has an amplitude between 50-100 mPa in the
medical range of dose delivery at 1-2 mGy/pulse. Based on the previous
experiments, the only conclusion we can make from this is that this amplitude
range is consistent with the upper limit of 0.2 Pa set by the experimental
amplitude of XA signal induced from a lead block Hickling et al. (2014b).
Nevertheless, an important notice in the analytically simulated therapeutic XA
signal as given in Fig 2 is that the relative amplitude between a crest and
trough is approximately the same for two boundaries consistently with the
experimental and k-wave simulated signals. Therefore, it can be further
normalized based on the experimental measurements of real pressure amplitudes
if necessary. Lastly, central frequency of the FFT spectrum of therapeutic XA
signal is found to be $f_{0}=45.8$ kHz consistently with the previous results
Xiang et al. (2013); Hickling et al. (2014b) with an error of 2.7%. Besides,
qualitative two-frequency content of the spectrum is similarly observed by
Xiang et al. (2013). The use of band-pass filter however filtered out the
higher frequency in Hickling et al. (2014b). An interesting but experimentally
unjustified observation based on our simulations is that the waveform
undergoes a phase shift and displacement as we move the transducer by
preserving relative pressure amplitude as in Fig 4(a) which might be accounted
by the dispersive nature of sound waves. These show that analytic approach can
be considered as an alternative to model therapeutic XA signals once the
explicit experimental verification which may bring further corrections is also
provided.

Since Bowen et al. (1991) observed the emission of induced acoustic signals
from water due to therapeutic X-ray dose gradient at the treatment edges,
succeeding XA studies Hickling et al. (2014a); Xiang et al. (2013); Hickling
et al. (2014b, 2016); Xiang et al. (2016) focused on imaging different objects
inside water-based mediums. Among these studies, only Xiang et al. (2013)
recorded therapeutic X-ray induced acoustic signal from the tissue by
operating 10 MV linac at 30 Gy/min which is a quite larger dose rate than
applied in radiotherapy, and Hickling et al. (2016) obtained therapeutic XA
signal from the water medium which can be also considered as a tissue model,
whereas other studies obtained XA signal from the irradiation of lead blocks
due to its larger density causing relatively large amplitude for the acoustic
signal. Physical magnitude of therapeutic XA pressure wave in tissue has still
been not experimentally established although Bowen et al. (1991) predicted
that it should be clearly visible in water at total doses less than 1 Gy. This
study provides for the first time with the quantitative parametrization of
radiotherapeutic XA waves. It verifies the feasibility of observing acoustic
signal from radiotherapeutic photon beams based on the pressure amplitudes and
NEP of an experimental transducer with an effective central frequency 50 kHz
and size $3\leavevmode\nobreak\ cm^{2}$.

This analytic exploration of pressure-absorption dependence can be useful both
for imaging and real-time dosimetry purposes, once XA signal is detected.
Putting X-ray imaging into clinical practice via XA waves requires the
minimization of the absorbed dose, and hence the X-ray beam energy is
restrained. Having the explicit relation between pulse duration and ultrasound
signal, a proper choice for the X-ray source can be made to obtain optimum
signal over the noise level. On the other hand, the fact that any fractional
change in absorption cross section of tissue compositions translates into
change in the induced pressure profile Xiang et al. (2016) can be used to
monitor whether the irradiated region is normal or cancerous. This might
enable simultaneous recognition of the amount of dose which is related with
the amplitude of the XA signal in healthy tissue and hence immediate
intervention to therapy session. Nevertheless, due to the uniform medium
assumption, absorption cross section is solely defined by the Gaussian width
parameter $\sigma$ determined by the field size on the surface by using the
fact that X-ray absorption coefficient is the inverse of the minimum dimension
of the area of heat evolution, in this study. This parameter can be important
to distinguish acoustic signals induced in mediums with different
characteristics. ‘Characteristics of the medium’/acoustic signal dependence is
previously pointed out by Xiang et al. (2016) although it has not been proven
by experiment or simulation. Depending on the medium characteristics,
probability of absorption is governed by the absorption cross section which is
an explicit parameter in our model. Hence, we argue that succeeding change of
the pressure profile due to the change in the width parameter $\sigma$ of the
Gaussian curve modeling the heating function $H(r,t)$ should mimic the
medium/XA signal dependence. Once this dependence is experimentally observed,
our claim, i.e. whether the dependence can be simulated by varying $\sigma$,
should be tested. On the other hand, we point out that any parameter
dependence other than the absorbed dose is implicit in the previous simulation
model Hickling et al. (2016), and hence established XACT simulation techniques
by no chance can simulate this dependence even if it will be experimentally
observed. Similar to the characterization of a bone with osteoporosis using
photoacoustic waveforms Feng et al. (2015), an experiment that measures XA
signals in a partially cancerous tissue might be considered for this purpose.

An analytic model that can fast and simply simulate acoustic waves induced by
radiotherapeutic X-rays is proposed for the first time. Considering the
previous studies on this subject, we believe that the significance of this
study is the foundation of a novel and self-contained analytic approach to
simulate the therapeutic X-ray acoustic waves based on the physical
parametrization of the energy transfer process. This not only provides a
better understanding of the physical phenomena underlying the medical
technique in terms of the medically relevant parameters such as field size,
pulse duration, absorbed dose per pulse etc. together with the physical
assumptions used to obtain a solution to the photo-acoustic equation, but also
brings consistent simulation results with previous experimental and k-Wave
results. This model can be easily implemented and modified for a quantitative
investigation of any X-ray acoustic study, and hence it is an alternative tool
to numerical approach. Further experimental research with more quantitative
analysis will improve physical modeling in this area.

## Disclosure of Conflicts of Interest

The authors have no relevant conflicts of interest to disclose.

## References

  * Bowen et al. (1991) T Bowen, CX Chen, SC Liew, WR Lutz, and RL Nasoni.  Observation of ultrasonic emission from edges of therapeutic x-ray beams.  _Physics in medicine and biology_ , 36(4):537, 1991.
  * Xiang et al. (2016) Liangzhong Xiang, Shanshan Tang, Moiz Ahmad, and Lei Xing.  High resolution x-ray-induced acoustic tomography.  _Scientific reports_ , 6, 2016.
  * Xiang et al. (2013) Liangzhong Xiang, Bin Han, Colin Carpenter, Guillem Pratx, Yu Kuang, and Lei Xing.  X-ray acoustic computed tomography with pulsed x-ray beam from a medical linear accelerator.  _Medical physics_ , 40(1), 2013.
  * Hickling et al. (2014a) S Hickling, M Hobson, and I El Naqa.  Feasibility of x-ray acoustic computed tomography as a tool for noninvasive volumetric in vivo dosimetry.  _International Journal of Radiation Oncology• Biology• Physics_ , 90(1):S843, 2014a.
  * Wang and Wu (2012) Lihong V Wang and Hsin-i Wu.  _Biomedical optics: principles and imaging_.  John Wiley & Sons, 2012.
  * Hickling et al. (2016) Susannah Hickling, Hao Lei, Maritza Hobson, Pierre Léger, Xueding Wang, and Issam El Naqa.  Experimental evaluation of x-ray acoustic computed tomography for radiotherapy dosimetry applications.  _Medical Physics_ , 2016.
  * Askaryan (1957) GA Askaryan.  Hydrodynamic radiation from the tracks of ionizing particles in stable liquids.  _The Soviet Journal of Atomic Energy_ , 3(2):921–923, 1957.
  * Hickling et al. (2014b) Susannah Hickling, Pierre Léger, and Issam El Naqa.  Simulation and experimental detection of radiation-induced acoustic waves from a radiotherapy linear accelerator.  In _Ultrasonics Symposium (IUS), 2014 IEEE International_ , 29–32. IEEE, 2014b.
  * (9) Lev D Landau and JB Sykes.  Fluid mechanics: Vol 6.  292–294.
  * Askariyan et al. (1979) GA Askariyan, BA Dolgoshein, AN Kalinovsky, and NV Mokhov.  Acoustic detection of high energy particle showers in water.  _Nuclear instruments and Methods_ , 164(2):267–278, 1979.
  * Bowen and Learned (1979) T Bowen and JG Learned.  Acoustic detection of ultra high energy neutrinos.  In _International Cosmic Ray Conference_ , volume 10, 386, 1979\.
  * Sulak et al. (1979) L Sulak, T Armstrong, H Baranger, M Bregman, M Levi, D Mael, J Strait, T Bowen, AE Pifer, PA Polakos, et al.  Experimental studies of the acoustic signature of proton beams traversing fluid media.  _Nuclear Instruments and Methods_ , 161(2):203–217, 1979.
  * Griffiths (2005) David J Griffiths.  Introduction to electrodynamics, 2005:423.
  * Learned (1979) John G Learned.  Acoustic radiation by charged atomic particles in liquids: An analysis.  _Physical Review D_ , 19(11):3293, 1979.
  * Baily (1992) Norman A Baily.  A review of the processes by which ultrasound is generated through the interaction of ionizing radiation and irradiated materials: Some possible applications.  _Medical physics_ , 19(3):525–532, 1992.
  * Jones et al. (2014) K C Jones, A Witztum, C M Sehgal, and S Avery.  Proton beam characterization by proton-induced acoustic emission: simulation studies.  _Physics in medicine and biology_ , 59(21):6549, 2014.
  * Karzmark and Pering (1973) CJ Karzmark and Neil C Pering.  Electron linear accelerators for radiation therapy: history, principles and contemporary developments.  _Physics in medicine and biology_ , 18(3):321, 1973.
  * Lyamshev (2004) Leonid M Lyamshev.  _Radiation acoustics_.  2004:41.
  * Wang et al. (2003) Xueding Wang, Yongjiang Pang, Geng Ku, Xueyi Xie, George Stoica, and Lihong V Wang.  Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain.  _Nature biotechnology_ , 21(7):803–806, 2003\.
  * Townson and Zavgorodni (2014) Reid W Townson and Sergei Zavgorodni.  Pre-treatment radiotherapy dose verification using monte carlo doselet modulation in a spherical phantom.  _Physics in medicine and biology_ , 59(8):1923, 2014.
  * Podgorsak (2005) Ervin B Podgorsak.  External photon beams: Physical aspects.  _Radiation Oncology Physics: A handbook for teachers and students’. Vienna: IAEA_ , 2005.
  * Ahmad et al. (2015) Moiz Ahmad, Liangzhong Xiang, Siavash Yousefi, and Lei Xing.  Theoretical detection threshold of the proton-acoustic range verification technique.  _Medical physics_ , 42(10):5735–5744, 2015.
  * Wang (2008) Lihong V Wang.  Tutorial on photoacoustic microscopy and computed tomography.  _IEEE Journal of Selected Topics in Quantum Electronics_ , 14(1):171–179, 2008.
  * Morse and Ingard (1968) Philip McCord Morse and K Uno Ingard.  _Theoretical acoustics_.  Princeton university press, 323-325, 1968.
  * Terunuma et al. (2007) Toshiyuki Terunuma, Takeji Sakae, Yoshinori Hayakawa, Akihiro Nohtomi, Yoshihisa Takada, Kiyoshi Yasuoka, and Akira Maruhashi.  Waveform simulation based on 3d dose distribution for acoustic wave generated by proton beam irradiation.  _Medical physics_ , 34(9):3642–3648, 2007.
  * Kruger et al. (1995) Robert A Kruger, Pingyu Liu, Yuncai Fang, C Robert Appledorn, et al.  Photoacoustic ultrasound (paus)—reconstruction tomography.  _Medical physics_ , 22(10):1605–1609, 1995.
  * Wang et al. (2004) Yi Wang, Da Xing, Yaguang Zeng, and Qun Chen.  Photoacoustic imaging with deconvolution algorithm.  _Physics in medicine and biology_ , 49(14):3117, 2004.
  * Maciszewski and Scharf (2004) Wieslaw Maciszewski and Waldemar Scharf.  Particle accelerators for radiotherapy: Present status and future.  In _Astroparticle, Particle and Space Physics, Detectors and Medical Physics Applications: Proceedings of the 8th Conference_ , page 402. World Scientific, 2004.
  * Mayles et al. (2007) Philip Mayles, Alan Nahum, and Jean-Claude Rosenwald.  _Handbook of radiotherapy physics: theory and practice_.  CRC Press, 2007.
  * Constantinou et al. (1982) Chris Constantinou, FH Attix, and Bhudatt R Paliwal.  A solid water phantom material for radiotherapy x-ray and $\gamma$-ray beam calibrations.  _Medical physics_ , 9(3):436–441, 1982.
  * Quach et al. (2000) KY Quach, J Morales, MJ Butson, Anatoly B Rosenfeld, and Peter E Metcalfe.  Measurement of radiotherapy x-ray skin dose on a chest wall phantom.  _Medical physics_ , 27(7):1676–1680, 2000.
  * Zhu et al. (1995) Yunping Zhu, Xun-Qing Jiang, and Jake Van Dyk.  Portal dosimetry using a liquid ion chamber matrix: Dose response studies.  _Medical physics_ , 22(7):1101–1106, 1995.
  * Butson et al. (1996) Martin J Butson, Anatoly Rozenfeld, Jagdish N Mathur, Martin Carolan, Tony PY Wong, and Peter E Metcalfe.  A new radiotherapy surface dose detector: the mosfet.  _Medical physics_ , 23(5):655–658, 1996.
  * O’Keeffe et al. (2016) Sinead O’Keeffe, Wenhui Zhao, Weimin Sun, Daxin Zhang, Zhuang Qin, Ziyin Chen, Yu Ma, and Elfed Lewis.  An optical fibre-based sensor for real-time monitoring of clinical linear accelerator radiotherapy delivery.  _IEEE Journal of Selected Topics in Quantum Electronics_ , 22(3):35–42, 2016.
  * Rogers et al. (1995) DWO Rogers, BA Faddegon, GX Ding, C-M Ma, J We, and TR Mackie.  Beam: a monte carlo code to simulate radiotherapy treatment units.  _Medical physics_ , 22(5):503–524, 1995.
  * Ding (2002) George X Ding.  Energy spectra, angular spread, fluence profiles and dose distributions of 6 and 18 mv photon beams: results of monte carlo simulations for a varian 2100ex accelerator.  _Physics in medicine and biology_ , 47(7):1025, 2002.
  * Jagetic and Newhauser (2015) Lydia J Jagetic and Wayne D Newhauser.  A simple and fast physics-based analytical method to calculate therapeutic and stray doses from external beam, megavoltage x-ray therapy.  _Physics in medicine and biology_ , 60(12):4753, 2015.
  * Taddei et al. (2013) Phillip J Taddei, Wassim Jalbout, Rebecca M Howell, Nabil Khater, Fady Geara, Kenneth Homann, and Wayne D Newhauser.  Analytical model for out-of-field dose in photon craniospinal irradiation.  _Physics in medicine and biology_ , 58(21):7463, 2013.
  * Winkler et al. (2013) Amy M Winkler, Konstantin Maslov, and Lihong V Wang.  Noise-equivalent sensitivity of photoacoustics.  _Journal of biomedical optics_ , 18(9):097003, 2013.
  * Koehler et al. (1975) AM Koehler, RJ Schneider, and JM Sisterson.  Range modulators for protons and heavy ions.  _Nuclear Instruments and Methods_ , 131(3):437–440, 1975.
  * Larsson (1961) Börje Larsson.  Pre-therapeutic physical experiments with high energy protons.  _The British journal of radiology_ , 34(399):143–151, 1961.
  * Clasie et al. (2010) Ben Clasie, Andrew Wroe, Hanne Kooy, Nicolas Depauw, Jay Flanz, Harald Paganetti, and Anatoly Rosenfeld.  Assessment of out-of-field absorbed dose and equivalent dose in proton fields.  _Medical physics_ , 37(1):311–321, 2010.
  * Bortfeld (1997) Thomas Bortfeld.  An analytical approximation of the bragg curve for therapeutic proton beams.  _Medical physics_ , 24(12):2024–2033, 1997.
  * Feng et al. (2015) Ting Feng, Joseph E Perosky, Kenneth M Kozloff, Guan Xu, Qian Cheng, Sidan Du, Jie Yuan, Cheri X Deng, and Xueding Wang.  Characterization of bone microstructure using photoacoustic spectrum analysis.  _Optics express_ , 23(19):25217–25224, 2015.


DRAFT

# A Multilayer Network Model implementation for Covid-19

Juan G. Calvo Universidad de Costa Rica, Centro de Investigación en Matemática
Pura y Aplicada - Escuela de Matemática, San José, Costa Rica
<EMAIL_ADDRESS><EMAIL_ADDRESS>luisalberto.barboza@ucr.ac.cr). Fabio Sanchez11footnotemark: 1 Luis A.
Barboza11footnotemark: 1 Yury E. García Universidad de Costa Rica, Centro de
Investigación en Matemática Pura y Aplicada, San José, Costa Rica
<EMAIL_ADDRESS>yury.garciapuerta@ucr.ac.cr). Paola
Vásquez22footnotemark: 2

###### Abstract

We present a numerical implementation for a multilayer network used to model
the transmission of Covid-19 or other diseases with a similar transmission
mechanism. The model incorporates different contact types between individuals
(household, social contacts, and strangers), which allows flexibility compared
to standard SIR type models. The algorithm described in this paper is a
simplification of the model used to give public health authorities an
additional tool for the decision-making process in Costa Rica, by simulating
extensive possible scenarios and projections.

###### keywords:

multilayer network, Covid-19, epidemic models, pandemic

###### AMS:

93A16, 92C60

##  1 Introduction

Throughout the Covid-19 pandemic, mathematical and statistical models have
been used as a way to forecast or analyze the behavior of the disease in
various parts of the world [1, 2, 3, 4, 5]. Different mathematical tools have
been implemented to understand the disease transmission and to plan for
strategic sanitary measures to minimize the burden on the public health system
[1, 5, 6]. Deterministic models have been widely used to model infectious
diseases in the past [7, 8, 9, 10, 11, 12]. However, due to human mobility,
unpredictable human behavior and other factors, adaptable models are better
suited to evolve with the ever-changing behavior of the population during the
pandemic as sanitary measures are implemented in a particular region of
interest.

In this paper, we consider an adaptive multilayer network, capable of managing
large populations efficiently; see [13, 14]. A probability of infection is
computed on a daily basis for each node independently, permitting a parallel
implementation for real-life applications. We count daily interactions per
individual and determine if the virus is transmitted depending on the
attributes of each interaction.

One of the main advantages of a network model is its flexibility when taking
into account individual and social behavior, that allows us to mimic sanitary
measures and evaluate possible impacts on the population. As the pandemic
progresses, it is possible to incorporate early or ongoing behavioral changes
such as reductions in community interactions, mask use, mobility restrictions,
network contacts, and lockdowns. We include a Matlab app with the algorithm
implementation that can be found in [15], where several attributes and
parameters can be modified111A full open code will be available as well..

The model we present in this paper is a simplification of the model that has
been developed and taken into consideration by the health authorities in Costa
Rica. It includes the entire population (approximately five million
individuals) and the 81222There are 82 municipalities in the country. However,
the 82nd was established in 2018 and there is no historical census data in
order to model mobility. municipalities of the country. The model also
considers different scenarios with gradual re-openings, mobility restrictions,
different percentages for mask use, as well as short and long-term projections
of the disease. Results were delivered to the public health authorities as
part of the decision-making process [16, 17, 18, 19]. The model takes into
account local data, public health policies pertinent to the country and more
specific details than the model described in this work.

The rest of the paper is organized as follows. In Section 2, we describe the
algorithm. A pseudo-code is presented in Section 3 with a discussion of the
implementation. In Section 4, we present two toy examples, in order to show
the behaviour of the model. Finally, we include some conclusions and future
work in Section 5.

##  2 Preliminaries

Consider a constant population that lives in a fixed number of counties.
Mathematically, we consider a graph (a set of nodes that are connected by
edges) where each node represents one individual. Interactions are classified
into three layers, which represent a different type of contact with particular
characteristics that account for the virus transmission. These layers are: (1)
a household network (people that live in the same house), (2) a social network
(known contacts such as friends and colleagues), and (3) a sporadic network
(strangers that you may encounter in short periods of time when you visit
random locations).

An edge between two nodes represents a feasible interaction between two
individuals at a particular time, with specific attributes for its interaction
depending on the layer. These layers are randomly generated. Layers 1 and 2
are fixed for each simulation and layer 3 can change for each time step, since
usually a node has no control on sporadic encounters.

For the first layer, it is assumed that we know the total number of households
per county, and the average number of individuals per household per county;
the information for Costa Rica is available at [20]. We then group all the
individuals into families. The number of members per household is chosen by
using a Poisson distribution with mean equal to the average number of
individuals per household per county. We assume total connectivity for each
family and no edges between different families. It is possible to randomly
assign an age group to each family member if data is available, in order to
include differences in interactions accordingly; see further details in
Section 3.1.

For the second and third layers, choosing the contacts of a node requires some
assumptions. It is natural to expect a different number of interactions
depending on age, location, density, socioeconomic factors and more. In our
implementation, for every node we need to: (1) define the degree of the node
(number of contacts) per layer, given by a uniform distribution on an interval
with endpoints equal to the minimum and maximum number of allowed contacts per
county per layer, and (2) choose randomly its contacts from different counties
based on a given connectivity distribution between counties; see Section 3.3.

Each node on the graph includes several attributes. We classify them in two
classes: (1) fixed parameters, such as county, household members, age group,
degree and graph connectivity for layer 2, and (2) variable attributes such as
connectivity for layer 3, epidemiological state, number of days at current
epidemiological state, number of interactions with social contacts per day,
mask use and self-care behaviour. The former are defined at the beginning of
each simulation, and the latter can take random values every time step
accordingly.

$S$$E$$O$$H$$R$$U$$D$ Fig. 1: Model compartment states and transitions.

For simplicity, the propagation in our model is based on an SIR type model
with seven compartment states defined by the following epidemiological
variables: susceptible ($S$), exposed ($E$), diagnosed or observed ($O$),
undiagnosed or not observed ($U$), hospitalized ($H$), recovered ($R$), and
dead ($D$) states; see Figure 1. Modifications with additional states and
transitions are straightforward to include.

SporadicSocialHousehold Fig. 2: An infectious node and its network.
SporadicSocialHousehold Fig. 3: Two new infected nodes.

SporadicSocialHousehold Fig. 4: Augmented network. SporadicSocialHousehold
Fig. 5: Four new infected nodes.

Fig. 6: Modeling the transmission of the virus. The network adapts on a daily
basis depending on new exposed nodes, imitating the evolution of new cases.

Once the multilayer network is defined, we model the transmission of the virus
by computing probabilities of infection $p_{i}$ for interactions between
infectious and susceptible nodes; see Section 3.4. For a given initial set of
exposed nodes, we consider their contacts via the different layers; see Figure
2. Depending on $p_{i}$, adjacent nodes of infectious nodes can become
exposed; see Figure 3. We then increase the network by including new exposed
nodes and their contacts as shown in Figure 4, continuing the process on a
daily basis; see Figure 5.

######  Remark 1.

The degree of each node in layer 2 is fixed for each simulation. Nevertheless,
on each time step a subset of edges is selected because only a small number of
daily contacts is expected to occur. As mentioned before, layer 3 changes
daily.

######  Remark 2.

If the population has a large number of individuals, it is not recommended to
create the whole graph initially due to execution times. Instead, at every
time step we add new nodes and their edges on the graph in order to achieve
faster running times; see Figure 6 and Section 3.3 for further details.

##  3 Implementation

In this section we present a pseudo-code for the modeling of Covid-19 as
described in Section 2; see Algorithm 1. We describe each step in the
following sections in a general and simplified setting.

Algorithm 1 Multilayer network model

1:procedure networkModel

2: Create household network. $\triangleright$ See Section 3.1

3: Define initial conditions.$\triangleright$ See Section 3.2

4: for day = $1,2,3\ldots$ do

5: for all new exposed nodes $E$ do $\triangleright$ See Section 3.3

6: if family has not been added then

7: Read connections in the family layer.

8: Store edge connectivity.

9: end if

10: if node requires more contacts in layer 2 then

11: Create list of eligible contacts.

12: Choose random nodes for layer 2.

13: Store edge connectivity for layer 2.

14: end if

15: end for

16: for all infectious nodes do

17: Create list of eligible contacts.

18: Choose random nodes for layer 3.

19: end for

20: for all susceptible nodes $S$ do $\triangleright$ See Section 3.4

21: Compute probability of infection.

22: Determine if current node has been infected.

23: end for

24: Determine transitions between states.$\triangleright$ See Section 3.5

25: Store results for current day.

26: end for

27:end procedure

###  3.1 Creation of household network

In this implementation we assume that the following data is available: total
number of households per county, average number of individuals per household
per county and age groups; see a toy example in Table 2. We then create all
the households and assign their inhabitants by using a Poisson distribution
with mean equal to the average number of individuals per household per county;
see Figure 7.

Each individual is labeled with an ID from $1$ to $N$, where $N$ is the size
of the population. We create the arrays IDtoCounty, IDtoFamily and IDtoAge
that return the county, family and age group of a given ID, and the cells
familyToIDs and countyToIDs that return the IDs for a given family or county,
respectively. These arrays allows us to directly access required information
to construct the graph.

Fig. 7: Three families with 4, 5 and 1 members. A fully connected graph for
each component is assumed. Colors represent different age groups.

###  3.2 Initial conditions

For simplicity, we assume an initial graph given by all the exposed
individuals at time zero and their contact networks; see Figure 2 for the case
of one initial exposed node. If data is available, we can include chains of
transmission for a fixed time as an initial condition, by creating edges among
identified clusters. Before community transmission in Costa Rica, the public
health authorities were able to keep track of chains of transmission for
approximately four months. This data allowed us to create accurate contact
networks for the early stages of the pandemic.

###  3.3 Updating the network

As mentioned in Remark 1, creating a particular graph and its structure for a
large population is both time and memory consuming. Instead, our graph evolves
over time. In order to keep track of the transmission of the disease, we only
need to include the edges between new exposed individuals $E$ and their
contacts.

We choose layer 2 contacts for new added nodes as follows. We first create a
list of eligible nodes based on a connectivity distribution between counties.
For instance, density, demographic information or mobility data could be used
for this purpose. For simplicity, we consider a binary symmetric matrix where
its entry $(m,n)$ is equal to $1$ if there are connections between counties
$m$ and $n$, and $0$ otherwise; see a toy example in Table 2. This means that
the eligible list for a node includes all IDs from counties that have
connectivity with its county. We then randomly choose the remaining number of
contacts for each node from the eligible list. Further restrictions can be
added when defining the eligible list, such as age group, density and distance
between counties. Edges on layer 3 are daily updated in a similar way,
excluding preassigned contacts.

###  3.4 Probability of infection

If the probability that node $j$ infects a susceptible node $i$ at a given day
is $\beta_{ij}$, then the probability of infection $p_{i}$ for node $i$ is
given by

$p_{i}=1-\prod_{j\neq i}(1-\beta_{ij})$

where $j$ includes all the indices of nodes that can infect node $i$ that
particular day. The value $\beta_{ij}$ can depend on the layer,
epidemiological state, social distancing, mask use and self-care behaviour,
among other factors. Transmission of the virus can also include imported cases
where infection occurs in an external location. Non-periodic events can also
be considered (as massive events or the creation of clusters due to super
spreader events) to study their impact and possible scenarios, where different
probabilities $p_{i}$ can be defined for each particular case.

In our implementation, we consider two binary random variables useMask and
selfCare that indicate whether or not an encounter between two nodes includes
mask use and self-care, respectively. We assume that a fixed fraction of
diagnosed cases do not isolate and that hospitalized individuals cannot infect
susceptible nodes. We then have different values for $\beta_{ij}$ depending on
the interaction between the two nodes and their attributes. We only need to
check susceptible nodes that have infectious contacts. Once $p_{i}$ is
determined, a uniformly distributed pseudo-random number $r_{i}\in(0,1)$ is
generated. If $r_{i}<p_{i}$, then node $n_{i}$ is marked as an exposed
individual at the current time step.

###  3.5 Transitions between states

At each time step, every node has one of the following states: susceptible,
exposed, diagnosed, undiagnosed, hospitalized, recovered or dead; see Figure
1. The event where a susceptible becomes exposed is given by the probability
of infection as discussed in Section 3.4. Any other event depends on a
probability of transition. Every node has a different probability to move to a
feasible next state depending on the period of time at its current
epidemiological state. If such event should happen, we then store the new
state and day; if not, the state remains unchanged.

We show a toy example of probabilities of transition for the diagnosed class
$O$ in Table 1. In this case, if a node has one day as $O$, it has a $40\%$
chance that it will be hospitalized. Otherwise, it will remain as $O$ on day
2. On day three, there are three possibilities: remain at $O$ ($50\%$), move
to $H$ ($10\%$) or recover $R$ ($40\%$). Note that each row adds up to 1 (one
state should always be chosen according to these probabilities), and
eventually all nodes that have been exposed should end up as $R$ or $D$. This
approach allows us to include variable stay periods on each class if data is
available; average stay periods can be considered otherwise.

time in $O$ | $S$ | $E$ | $O$ | $U$ | $H$ | $R$ | $D$
---|---|---|---|---|---|---|---
1 |  |  | 0.6 |  | 0.4 |  |
2 |  |  | 1.0 |  |  |  |
3 |  |  | 0.5 |  | 0.1 | 0.4 |
4 |  |  | 0.1 |  | 0.3 | 0.6 |
5 |  |  |  |  |  | 1.0 |
Table 1: Toy example for the probability of transitions between states from
the observed class $O$; omitted entries are zero. From compartment $O$, it is
possible to move to states $H$ or $R$ as shown in Figure 1.

##  4 Numerical results

In this section we present results for two toy examples obtained with a
parallel code implemented in Matlab and a Lenovo SR650 Server with two 2.20GHz
Intel Xeon Plata 4214 processors. We remark that a generalization of this
implementation has allowed us to study disease projections and different
scenarios in Costa Rica, where the total population is approximately five
million people with feasible running times.

###  4.1 Fixed parameters

We consider four counties with a population of 1300 nodes; its demographic
distribution and the connectivity matrix are shown in Table 2. We consider
initially $E_{0}$ exposed nodes randomly chosen. Parameters used in these
simulations appear in Table 3. We assume that $70\%$ of the population wears a
face mask, $35\%$ respects self-care and $60\%$ of exposed nodes become
diagnosed. Values for the probability of infection are shown in Table 3 for
layer 2; probabilities are halved for layer 3. Finally, the probabilities of
transitions between states are shown in Appendix 1.

County | HH | ANIHH | AG1 | AG2 | AG3
---|---|---|---|---|---
1 | 140 | 3.0 | 20 | 50 | 30
2 | 100 | 3.5 | 28 | 44 | 28
3 | 80 | 4.0 | 20 | 40 | 40
4 | 60 | 3.5 | 40 | 40 | 20

County | 1 | 2 | 3 | 4
---|---|---|---|---
1 | 1 | 1 | 1 | 1
2 | 1 | 1 | 0 | 1
3 | 1 | 0 | 1 | 0
4 | 1 | 1 | 0 | 1

Table 2: (left) Toy demographic distribution considered in Section 4.1. We
present the number of households (HH), the average number of individuals per
household (ANIHH) and age groups (AG) percentage distribution. (right)
Connectivity matrix. Non-zero entries indicate possible connections among
counties.

Parameter | Interval
---|---
Contacts in social network | $[5,15]$
Contacts in sporadic network | $[0,15]$
Contacts per day, county 1 | $[5,15]$
Contacts per day, county 2 | $[3,13]$
Contacts per day, county 3 | $[1,10]$
Contacts per day, county 4 | $[1,10]$

Face mask | Self care | $\beta$
---|---|---
No | No | 0.21
No | Yes | 0.15
Yes | No | 0.08
Yes | Yes | 0.05

Table 3: (left) Network parameters in toy example and (right) probability of
infection depending on the type of encounter between two nodes for layer 2.
Values on layer 3 are halved; see Section 4.1.

We present a simulation for $E_{0}=1$ (only one initial exposed node) in
Figure 8. In Figure 9, we present 1000 curves for the cumulative number of
cases. We remark that in this scenario with only one initial exposed node,
there is no outbreak in $7\%$ of the simulations. Meanwhile, in the remaining
$93\%$ there are 700 cumulative number of cases on average.

Fig. 8: One simulation

Fig. 9: Cumulative number of cases

Fig. 10: Experiments with $E_{0}=1$. (left) Cumulative number of cases, $O$,
$U$ and $H$ for one simulation. (right) Cumulative number of cases for 1000
simulations; see Section 4.1. The variability is due to the fact that there is
only one initial exposed node.

For $E_{0}=10$ and the given parameters, there is a one-wave pattern in the
numerical results. We present the mean of 1000 simulations with their $5^{\rm
th}$ and $95^{\rm th}$ percentiles in Figure 15. Results can also be
desegregated by county and/or age group, but we omit these figures for the
sake of brevity. They can be obtained by using the app given in [15].

Fig. 11: Cumulative number of cases.

Fig. 12: Observed individuals.

Fig. 13: Hospitalized individuals.

Fig. 14: Undiagnosed individuals.

Fig. 15: Mean (black line) for 1000 simulations with $E_{0}=10$. The shaded
region corresponds to values between the $5^{\rm th}$ and $95^{\rm th}$
percentiles of 1000 simulations; see Section 4.1.

###  4.2 Modeling changes in behaviour

We now consider a population of one million individuals distributed in 15
counties. We impose restrictions on social behaviour by increasing the
percentage of individuals that wear a mask and decreasing the number of daily
interactions on an period of two months, starting on day 20, emulating the
effect of mild regulations. We present results for 1000 simulations in Figure
20, where we compare the effect of such restrictions.

Fig. 16: Cumulative number of cases.

Fig. 17: Observed individuals.

Fig. 18: Hospitalized individuals.

Fig. 19: Undiagnosed individuals.

Fig. 20: Percentage of cases with no measures (dashed line) and mild
restrictions (solid line) for two months starting on day 20. We consider a
reduction on the number of daily contacts and an increment in use of masks
with $E_{0}=10$. The shaded region corresponds to values between the $5^{\rm
th}$ and $95^{\rm th}$ percentiles of 1000 simulations; see Section 4.2.

##  5 Conclusions and final remarks

We have presented an implementation for a multilayer network to support public
health authorities during the Covid-19 pandemic. Capturing the transmission
dynamics of the SARS-CoV-2 virus is essential in order to provide information
about disease projections and specific data-driven scenarios.

The model exhibits a flexible structure that allowed us to include particular
attributes for each node and their contact networks, describing the
complexities of social behavior and the specific disease. In this article, the
implementation of the general structure of the model is described which can be
used in other social and epidemic contexts.

With a health emergency such as the one caused by the SARS-CoV-2 virus, public
health authorities were forced to make short-term decisions. Mathematical
models help authorities to systematize the transmission mechanisms of the
virus and project the expected behavior in the short and medium-term under
some assumptions. In particular, a network model is a valuable tool as an
input in the decision-making process for health authorities. It has the
capability to include social factors and an epidemiological structure,
enabling a greater level of detail in the study of the disease. Despite the
flexibility of the model, there are limitations that mostly pertain to data
availability and its quality.

Implementing and executing these models can be challenging. The availability
of information is closely related to the specificity of the question that
needs to be answered. Information on the initial conditions for the model and
parameter value estimates can be a limiting factor when data is not available
or its quality is not optimal. Efficient and fast algorithms are necessary due
to the on-demand health scenarios that need to be simulated frequently, and
execution times are limited by computational resources.

In an rapidly evolving epidemic as Covid-19, the model must be calibrated
frequently, as social behavior is constantly changing and, in many cases, it
is unpredictable. The model we have discussed in this paper allows us to
straightforwardly adapt and include frequent changes on variables that mimic
social behavior. Additionally, the model can be generalized to include other
factors, such as the effect of herd immunity and the role of vaccination, as
countries take the necessary steps to get back to the new normal.

## Acknowledgments

The authors would like to thank support by the Research Center in Pure and
Applied Mathematics and the Department of Mathematics at Universidad de Costa
Rica.

## References

  * [1] Center for Disease Control. COVID-19 Forecasts: Deaths — CDC. https://www.cdc.gov/coronavirus/2019-ncov/covid-data/forecasting-us.html, February 2021. (Accessed on 02/17/2021).
  * [2] Adiga, A., Dubhashi, D., Lewis, B., Marathe, M., Venkatramanan, S., Vullikanti, A. Mathematical models for Covid-19 pandemic: a comparative analysis. _Journal of the Indian Institute of Science_ , 100:793-807, 2020.
  * [3] Torrealba-Rodriguez, O., Conde-Gutiérrez, R.A., Hernández-Javier, A.L. Modeling and prediction of COVID-19 in Mexico applying mathematical and computational models. _Chaos, Solitons & Fractals_, 138:109946, 2020.
  * [4] Kucharski, A.J., Russell, T.W., Diamond, C., Liu, Y., Edmunds, J., Funk, S., Eggo, R.M., Sun, F., Jit, M., Munday, J.D. et al. Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The lancet infectious diseases, 20(5):553-558, 2020.
  * [5] Centre for the Mathematical Modelling of Infectious Diseases LSHTM. https://www.lshtm.ac.uk/research/centres/centre-mathematical-modelling-infectious-diseases, February 2021. (Accessed on 02/17/2021).
  * [6] Star, L., Moghadas, S.M. The role of mathematical modelling in public health planning and decision making. _Purple Paper, National Collaborative Center for Infectious Diseases_ , 2010.
  * [7] Perelson, A.S., Ribeiro, R.M. Modeling the within-host dynamics of HIV infection. _BMC biology_ , 11(1):1-10, 2013.
  * [8] Dietz, K., Schenzle, D. Mathematical Models for Infectious Disease Statistics. In: Atkinson A.C., Fienberg S.E. (eds), A Celebration of Statistics. Springer, New York, NY. https://doi.org/10.1007/978-1-4613-8560-8_8.
  * [9] Etbaigha, F., Willms, A.R., Poljak, Z. An SEIR model of influenza A virus infection and reinfection within a farrow-to-finish swine farm. _PLoS One_ , 13(9):e0202493, 2018.
  * [10] Hethcote, H.W. The mathematics of infectious diseases. _SIAM Review_ , 42(4):599-653, 2000.
  * [11] Choisy, M., Guégan, J.-F., Rohani, P. Mathematical modeling of infectious diseases dynamics. _Encyclopedia of infectious diseases: modern methodologies_ , 379, 2007.
  * [12] Chen, S.-C., Hsieh, M.-H. Modeling the transmission dynamics of dengue fever: implications of temperature effects. _Science of the total environment_ , 431:385-391, 2012.
  * [13] Kivelä, M., Arenas, A., Barthelemy, M., Gleeson, J.P., Moreno, Y., Porter, M.A. Multilayer networks. _Journal of Complex Networks_ , 2(3):203-271, 2014.
  * [14] De Domenico, M., Granell, C., Porter, M.A., Arenas, A. The physics of spreading processes in multilayer networks. _Nature Physics_ , 12(10):901-906, 2016.
  * [15] EpiMEC. Covid-19 app. www.github.com/epimec, 2021.
  * [16] Pan American Health Organization (PAHO/WHO). Costa Rica: Pandemia covid-19. Informe estratégico mensual $N^{o}$3\. https://www.paho.org/es/documentos/costa-rica-pandemia-covid-19-informe-estrategico-mensual-no-3, 2020. (Accessed on 03/09/2021).
  * [17] Pan American Health Organization (PAHO/WHO). Costa Rica: Pandemia covid-19. Informe estratégico mensual $N^{o}$4\. https://www.paho.org/es/documentos/costa-rica-pandemia-covid-19-informe-estrategico-mensual-no-4, 2020. (Accessed on 03/09/2021).
  * [18] Pan American Health Organization (PAHO/WHO). Costa Rica: Pandemia covid-19. Informe estratégico mensual $N^{o}$5\. https://www.paho.org/es/documentos/costa-rica-pandemia-covid-19-informe-estrategico-mensual-no-5, 2020. (Accessed on 03/09/2021).
  * [19] Pan American Health Organization (PAHO/WHO). Costa Rica: Pandemia covid-19. Informe estratégico mensual $N^{o}$6\. https://www.paho.org/es/documentos/costa-rica-pandemia-covid-19-informe-estrategico-mensual-no-6, 2020. (Accessed on 03/09/2021).
  * [20] INEC. Población — Instituto Nacional de Estadística y Censos. https://www.inec.cr/poblacion. (Accessed on 02/17/2021).

##  Appendix A Probabilities of transition

We present the probabilities of transition between epidemiological states used
in the toy examples of Section 4. We remark that there are different
probabilities that depend on the age group; see Tables 4, 5, 6, 7.

Days | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
$E$ | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | 0.8 |
$O$ |  |  |  |  |  | 0.1 | 0.6
$U$ |  |  |  |  |  | 0.1 | 0.4
Table 4: Probabilities from $E$ to $O$ and $U$ depending on the number of days at state $E$; empty entries are zero. Days | 1 | … | 12 | 13 | 14 | 15 | 16 | 17 | 18
---|---|---|---|---|---|---|---|---|---
$U$ | 1.0 | … | 1.0 | 0.9 | 0.8 | 0.7 | 0.5 | 0.1 |
$R$ |  |  |  | 0.1 | 0.2 | 0.3 | 0.5 | 0.9 | 1.0
Table 5: Probabilities from $U$ to $R$ depending on the number of days at state $E$; empty entries are zero. | Age group 1 | Age group 2 | Age group 3
---|---|---|---
Days | $O$ | $H$ | $R$ | $O$ | $H$ | $R$ | $O$ | $H$ | $R$
1 | 1.000 |  |  | 1.000 |  |  | 1.000 |  |
2 | 1.000 |  |  | 1.000 |  |  | 1.000 |  |
3 | 1.000 |  |  | 1.000 |  |  | 1.000 |  |
4 | 0.978 | 0.022 |  | 0.953 | 0.047 |  | 0.629 | 0.371 |
5 | 1.000 |  |  | 1.000 |  |  | 1.000 |  |
⋮ | ⋮ |  |  | ⋮ |  |  | ⋮ |  |
17 | 1.000 |  |  | 1.000 |  |  | 1.000 |  |
18 |  |  | 1.000 |  |  | 1.000 |  |  | 1.000
Table 6: Probabilities from $O$ to $H$ and $R$ depending on the number of days at state $E$; empty entries are zero. | Age group 1 | Age group 2 | Age group 3
---|---|---|---
Days | $H$ | $R$ | $D$ | $H$ | $R$ | $D$ | $H$ | $R$ | $D$
1 | 1.00 |  |  | 1.00 |  |  | 1.00 |  |
2 | 1.00 |  |  | 1.00 |  |  | 1.00 |  |
3 | 1.00 |  |  | 1.00 |  |  | 1.00 |  |
4 | 1.00 |  |  | 1.00 |  |  | 1.00 |  |
5 | 1.00 |  |  | 1.00 |  |  | 1.00 |  |
6 | 1.00 |  |  | 1.00 |  |  | 0.75 |  | 0.25
7 | 1.00 |  |  | 1.00 |  |  | 0.75 |  | 0.25
8 | 1.00 |  |  | 1.00 |  |  | 0.75 |  | 0.25
9 | 1.00 |  |  | 1.00 |  |  | 1.00 |  |
10 |  | 1.00 |  | 1.00 |  |  | 1.00 |  |
11 |  | 1.00 |  |  | 0.96 | 0.04 | 1.00 |  |
12 |  | 1.00 |  |  | 0.96 | 0.04 | 1.00 |  |
13 |  | 1.00 |  |  | 0.96 | 0.04 | 1.00 |  |
14 |  | 1.00 |  |  | 0.96 | 0.04 |  | 0.75 | 0.25
Table 7: Probabilities from $H$ to $R$ and $D$ depending on the number of days
at state $H$; empty entries are zero.