# An Insight of Heart-Like Systems with Percolation

Md Aquib Molla  0000-0003-0416-1349 Vidyasagar College, 39 Sankar Ghosh lane,
Kolkata 700006, India.  Sanchari Goswami  0000-0002-4222-5123 Vidyasagar
College, 39 Sankar Ghosh lane, Kolkata 700006, India.

###### Abstract

We study the signal percolation through heart-like biological system. Starting
from an initial distribution of Waiting and Inactive cells with probabilities
$p$ and $(1-p)$ respectively the signal propagation is observed in terms of
Active cells. As the signal enters the system from one end, the number of
arrival of Active sites at the other end is studied and analysis of the system
behaviour is made by varying a few important parameters of the system like
$p_{switch}$ (switching probability from Inactive to Waiting) and $p_{act}$
(switching probability from Waiting to active). In this connection, the non-
regular heart rhythms are discussed. Fraction of percolating paths through the
system is studied and showed a transition from $0$ to $1$ near $p=p_{c}$. Some
other important quantities like tortuosity and cluster distribution are
discussed and corresponding exponents have been found.

## Introduction

Percolation is one of the simplest and fundamental models in statistical
mechanics which is not always an exactly solvable model. In spite of its
simple update rules, it can be applied to describe a wide variety of systems.
Not only for the study of the dynamics of fluids through porous media [1] but
it can as well be applied to study the electrical movement in neurons,
fibrosis in organs, damage spreading in solids through fiber bundle model,
disease spreading in a community, opinion formation in society or even for
knowledge percolation [2, 3, 4, 5, 6, 7, 8, 9].

A very popular model of Percolation is the forest fire model (FFM) introduced
in [10, 11] mainly as a 2D lattice. There are several variations of the model
[11, 12]. In one, the fire can spread from top to bottom and from left to
right can be regarded as a general “Directed Percolation” (DP) model. For DP,
it is well understood that the bottom to top and right to left movements are
prohibited [13, 14]. There are other models of percolation where all top to
bottom, bottom to top, left to right and right to left movements are all
allowed and can be termed as Isotropic Percolation (IP). However, the model
where the percolation is allowed in the top to bottom direction and sideways
can be termed as Semi Directed Percolation (SDP). This was introduced for the
first time [15] in a random resistor-diode network in [16]. It has also been
studied in square and triangular lattices [17]. However in both [16, 17], it
was termed not as SDP but as Partially Directed Percolation.

The movement of signal in neurons for the heart and brain systems [2, 3, 4, 5]
can be studied in the light of percolation. Movement of heat and electrical
signal through different media have been studied in [19, 18] by percolation
methods. In addition to this, there are certain important quantities like
‘Tortuosity’ which have been studied for a few systems in combination with
percolation [20].

This process of signal percolation through heart is an interesting topic to
study. When the heartbeats do not follow the usual pattern but become very
fast or irregular, it leads to a disease called arrhythmia [21]. Arrhythmias
in turn cause health situations such as fainting, stroke, heart attack or even
death. Normally, a special group of cells in the sinoatrial (SA) node in the
upper right chamber of the heart begin the electrical signal to start the
heartbeat. This signal travels down the heart towards the ventricles, i.e.,
the two lower chambers of the heart. This organized pattern helps the heart to
beat in its usual way. However, several problems affect the normal conduction
of signal and lead to abnormal heart rhythms. Due to faster beats, there may
not be enough time in between beats and this may prevent the heart from
pumping the required amount of blood to the whole body which may lead to
arrhythmia and related health problems.

Several biological models have been used to study the heart system with the
help of computation. Computational models in cardiology showed enormous scope
in revealing diagnostic information related to cardiac disease. At the same
time there is provisions that these models help the diagnosis to be cost-
effective and to have low risk [5]. The cells can be activated by an external
stimulus which initiates an Action Potential (AP). The cells after activation
will go through a certain refractory period. The exact physiological processes
are described in [22] in a detailed manner. In [6], the model of heart is
studied with three states of cells: active, inactive and waiting in connection
with the AP.

One of the simplified system to study a percolation problem is a $L\times L$
lattice, where there is a parameter called percolation probability $p$ with
which a site could be open, i.e. can percolate “something”. This in turn
indicates that a site is closed with $1-p$ probability, means unable to
percolate. It is well known from different studies that below a critical
percolation probability $p_{c}$ the system is unable to percolate, whereas
above $p_{c}$ infinite percolation cluster can be formed. In connection to
percolation there are a number of exponents which can be studied as well.
Those exponents actually describes the system near criticality.

## Description of the Model

We study a 2D square grid of $L\times L$ cells where initially the cells may
be in two states, “Waiting” and “Inactive” as described in [6] also. A
“Waiting” cell is an open cell which is ready to hold anything percolating
through the system. An “Inactive” cell is however equivalent to a closed cell
which is unable to percolate. We choose the system with an initial
distribution of “Waiting” cells with probability $p$ and “Inactive” cells with
probability $1-p$. A “Waiting” cell can be activated with an action potential
(the thing that is actually percolating through the system) and then it
becomes an “Active” one. The condition for a Waiting cell to transform into an
Active one in the next step is that one or more of its “nearest cell
neighbors” should be Active. An “Active” cell which is holding the Action
Potential at a particular instant will be transformed into an “Inactive” one
in the next time step. For a heart system this is consistent with the fact
that there is a certain refractory period of the cells.

From now on, we will use numbers to indicate the type of the cells for our
convenience. Thus

  * •

an “Waiting” cell will be identified as “$0$”,

  * •

an “Active” cell will be identified as “$1$”,

  * •

an “Inactive” cell will be identified as “$2$”.

The general rule of transformation for a cell may be described as: Waiting
$\rightarrow$ Active $\rightarrow$ Inactive, or, using the numbers as
$0\rightarrow 1\rightarrow 2$. However, in some specific cases an “Inactive”
cell may become a “Waiting” one, i.e., $2\rightarrow 0$ is also allowed. This
is similar to the tree regrowth at several sites in Forest Fire Model as
discussed in [10, 11].

In addition to these rules, some definite probabilities are considered to be
associated with this system. We introduce two probabilities:

  * •

$p_{act}$(inhibitory) by which an Waiting cell can become an Active one, i.e.,
$0\rightarrow 1$.

  * •

$p_{switch}$(refractory) by which an Inactive cell can be reverted to a
Waiting one, i.e., $2\rightarrow 0$.

.

The AP is initiated from one side (here from the top) of the 2D matrix. This
is done by labelling all the cells on the first row as $1$. The different rows
are updated at different times $t$. It is obvious that to receive the AP at
the bottom/final row the number of time steps required will at least be equal
to the system size $L$.

The system is studied from two aspects. First, the number of arrival $N_{A}$
at the other end of the 2D matrix is observed with variation of the parameters
$p_{act}$ and $p_{switch}$ in subsection A of the Results. If there is no
Inactive cell, the AP reach all the cells of the final row in time $L$.
However, when there is a distribution of $0$ and $2$ cells, all the cells of
the final row do not contain the AP at the same time but $1$s arrive at
different cells at different times. At the final row the number of cells
containing the AP is counted at each time. Second, the percolation model of
the system is analyzed. The delay in arrival is measured and some interesting
properties related to that are studied. In the same context the cluster
distribution and several exponents related to the system are also studied. The
corresponding results are shown in subsections B and C of the Results.

## Results

###  .1 Study of Number of Arrival

The Action Potential as initiated from one end of the 2D lattice of size $L$
percolates through the system. The number of arrivals at the other end, i.e.,
the active sites at that ends $N_{A}$ is checked with time $t$ . We will study
the number of arrivals $N_{A}$ with two conditions:

  1. 1.

When the $2\rightarrow 0$ transition is not allowed $p_{switch}=0.0$ : This
arises when the cells after becoming Active receiving the AP for once become
inactive and stay in their refractory period for a very long time. Here
$p_{switch}=0.0$ is equivalent to an infinite refractory period. This can be
often seen due to certain health condition. For old people this is very common
[23]. At the same time it may also happen that the transition from
$0\rightarrow 1$ can occur with probability less than $1$. This may also
happen due to some irregularity in the cell activation [24].

  2. 2.

When $2\rightarrow 0$ transition is allowed with a certain probability other
than $0.0$ : Here $p_{switch}\neq 0.0$ indicates that the cells are in their
refractory period. However, very small value of $p_{switch}$ indicates that
the cells are in their refractory period for longer than usual. Medications
and better lifestyle may lead the system to have a higher $p_{switch}$.
However, having a very larger value of $p_{switch}$ is not suitable for this
case and will be discussed later.

In addition to $p_{switch}=0.0$ there may be variation of $p_{act}$ with which
we have studied the system.

The variation of number of arrival $N_{A}$ with $t$ is shown in Fig. 1 for
$p_{act}=1$ and $p_{switch}$ varying from $0.0$ to $0.9$. The plots show some
kind of oscillatory behaviour for high values of $p_{switch}$ and finally
saturate to some specific non-zero values $(N_{A})_{sat}$.

Figure 1: $N_{A}$ vs $t$ with $L$ = 200, $p=0.70$, $p_{act}=1.0$ and
$p_{switch}$ varying from $0$ to $0.90$. In the inset $p_{switch}=0.90$ has
been shown with equation 3 where $(N_{A})_{sat}\simeq 54$, $A=55$ and
$B=0.05$.

The oscillatory behaviour can be explained in the following way. It is well
known that when a fluid is flowing through a porous medium, there is some kind
of relative motion [25] which generates a fluid shear [19], can be considered
as a friction/damping. In the same way, when an electrical signal, here the
AP, is moving through the heart, there must be some retardation present. For
high $p_{switch}$, the motion is similar to an underdamped oscillatory one.
This is obvious as for high $p_{switch}$, the signal more frequently finds a
path to percolate, be it in forward or backward directions or sideways.
Therefore number of arrival $N_{A}$ is never zero. However, for
$p_{switch}=0.0$, it has been observed that the behaviour of $N_{A}$ is
Gaussian. This signifies that for low $p_{switch}$ the number of active cells
die down almost completely in time $L<t<2L$. Therefore for further AP signals
there is no chance to pass down through the system. However, for higher
$p_{switch}$, this is not true and even after a very long time there is always
a chance for the next APs to pass through the system. This is because here we
have $(N_{A})_{sat}>0$ always that indicateds a certain finite number of
Active cells will always be present in the system. Even if we study the system
for $t\rightarrow\infty$ the system still shows non-zero active sites.

It is already observed that for $p_{switch}\rightarrow 1$ the system shows a
damped oscillatory behaviour. Therefore it is worth comparing the system
parameters to the parameters associated with a damped oscillator. The
differential equation of a damped oscillator is

$\frac{d^{2}x}{dt^{2}}+2B\frac{dx}{dt}+\omega_{0}^{2}x=0$ (1)

where the coefficient $2B$ is related to damping and $\omega_{0}$ is the
natural frequency of oscillation. For underdamped/oscillatory motion
$B^{2}<\omega_{0}^{2}$ and the corresponding solution can be written as

$x(t)=A\cos(\omega t)\exp(-Bt)$ (2)

where $\omega=\sqrt{\omega_{0}^{2}-B^{2}}$. For this purpose, we fit $N_{A}$
with the following function:

$N_{A}(t)=(N_{A})_{sat}+A\cos(\omega{t^{\prime}}^{\nu})\exp(-Bt^{\prime})$ (3)

where $t^{\prime}=t-L$. The approximate values are shown in Table 1.

$p_{switch}$ |  $(N_{A})_{sat}$ |  $A$ |  $\omega$ |  $\nu$ |  $B$
---|---|---|---|---|---
0.90 | 53.83 | 55 | 1.00 | 1.10 | 0.05
0.80 | 48.78 | 58 | 0.65 | 1.20 | 0.07
0.70 | 45.31 | 70 | 0.45 | 1.29 | 0.09
0.60 | 41.21 | 75 | 0.31 | 1.37 | 0.10
Table 1: Approximate values of the parameters in Equation 3 with $L=200$.

In case of $p_{act}=p_{switch}=1$ (not shown in Fig. 1) the oscillatory
behaviour continues as the rows are updating as waiting $\rightarrow$ active
$\rightarrow$ inactive $\rightarrow$ waiting.

We now fix $p_{switch}=0$ and study the system behaviour by varying $p_{act}$.
It is already clear from Fig. 1 that the behaviour will be approximately
Gaussian. In Fig. 2 the behaviour is shown for $0.85\leqslant p_{act}\leqslant
1.0$. The parameters related to the Gaussian behaviour is shown in Table 2 for
$L=200$ by fitting with the following function:

$N_{A}(t)=(N_{A})_{max}e^{-\frac{1}{2}\left(\frac{t-\mu}{\sigma_{sd}}\right)^{2}}.$
(4)

where $(N_{A})_{max}$ is the peak value, $\mu$ is the mean and $\sigma_{sd}$
is the standard deviation of the Gaussian profile.

Figure 2: $N_{A}$ vs $t$ with $L=200$, $p_{switch}=1.0$ and $p_{act}$ varies from $0.85$ to $1.00$. $p_{act}$ |  $(N_{A})_{max}$ |  $\mu$ |  $\sigma_{sd}$
---|---|---|---
$1.00$ | $7.78$ | $259.10$ | $6.35$
$0.97$ | $6.20$ | $268.32$ | $7.48$
$0.95$ | $4.89$ | $275.54$ | $9.22$
$0.90$ | $2.66$ | $303.00$ | $14.01$
$0.85$ | $0.66$ | $356.50$ | $33.71$
Table 2: Approximate values of the parameters in Equation 4 with $L=200$.

It is observed that as $p_{act}$ decreases, $(N_{A})_{max}$ decreases. It is
in agreement to our description of the model. For $p_{act}<1$ the Waiting
cells do not receive the AP with $100\%$ probability and therefore cannot
become Active. Therefore $(N_{A})_{max}$ decreases. The movement of the peak
towards right with decreasing $p_{act}$, i.e., shifting of the mean of the
Gaussian profile towards higher $t$ can also be explained in the following
manner: As $p_{act}<1$ the signal cannot find its way smoothly inside the
system and therefore must take a longer time to appear at the other end. The
variations of $(N_{A})_{max}$, $\mu$ and $\sigma_{sd}$ are shown as a function
of $p_{act}$ in Fig, 3 for $L=50,100,200$.

Figure 3: Variations of the Gaussian parameters with $p_{act}$ where
$L=50,100$ and $200$, (a) $\frac{(N_{A})_{max}}{L^{\eta_{A}}}$ versus
$p_{act}$ with $\eta_{A}=0.7$, $m=1.16\pm 0.05$ and $c=-0.98\pm 0.05$, (b)
$\frac{\mu}{L^{\eta_{\mu}}}$ versus $p_{act}$ with $\eta_{\mu}=1.0$,
$a_{\mu}=1.26\pm 0.03$ and $b_{\mu}=-2.02\pm 0.20$ and (c)
$\frac{\sigma}{L^{\eta_{\sigma}}}$ vs $p_{act}$ with $\eta_{\sigma}=1.0$,
$a_{\sigma}=0.96\pm 0.03$ and $b_{\mu}=-7.82\pm 0.36$.

From Fig. 2 it is evident that with decrease of $p_{act}$ the curve broadens
and the height of the peak decreases. This suggests that not only arrivals are
delayed because of obstacles within the system, the signal is getting weakened
too. If we analyze the area under the curve carefully, it is evident that it
is decreasing with decrease in $p_{act}$. This also confirms presence of some
kind of retardation as it usually present for a real heart system.

###  .2 Study of Percolating Paths and Tortuosity

Next we study the probability $P_{N}$ in terms of the fraction of percolating
paths out of a $N_{T}$ number of trials against $p$. This is shown in Fig. 4
for $p_{act}=1.0,p_{switch}=0.0$ for different $L$. Here we have taken
$N_{T}=1000$. In the inset of Fig. 4 the data collapse has been shown. It has
been observed that the percolation threshold is at $p_{c}\simeq 0.59$ which is
in good agreement with [26, 27, 28, 29]. The function is checked to be fitted
with the function $\frac{(1+\tanh{\kappa x})}{2}$ where $\kappa=1.4$ and
$x=\frac{p-p_{c}}{p_{c}}$. The same has also been studied for other
combinations of $p_{act}$ and $p_{switch}$. As we decrease the value of
$p_{act}$, the percolation threshold $p_{c}$ shifts towards higher side as the
Waiting sites are becoming Active with less probability. If keeping
$p_{act}=1.0$ we increase the value of $p_{switch}$ even slightly, $p_{c}$
shifts towards lower $p$. This is also in agreement with Fig. 1. These are
shown in Fig.5 (a) for $p_{act}=0.80$, $p_{switch}=0.00$ and (b)
$p_{act}=1.0$, $p_{switch}=0.03$. From Fig. 5(a) $p_{c}$ is found to be $0.72$
and from 5(b) it is found to be close to $0.12$. It is to be noted that in
Fig. 5(b) there is no crossing and even for very small $p_{switch}$, $p_{c}$
drops rapidly.

Figure 4: $P_{N}$ vs $p$ with $p_{act}=1.0$, $p_{switch}=0.00$ where
$L=20,30,50,70$ and $100$ in the inset data collapse has been shown as
$\frac{(1+\tanh(\kappa x))}{2}$ with $n=1.5$, and $\kappa=1.4$ Figure 5:
$P_{N}$ vs $p$ where $L=20,30,50,70$ and $100$ with $n=1.5$, $(a)$
$p_{act}=0.80$ and $p_{switch}=0.00$ in the inset data colaps has been shown
as $\frac{(1+\tanh(\kappa x))}{2}$, $n=1.5$ and $\kappa=1.5$ $(b)$
$p_{act}=1.0$ and $p_{switch}=0.03$ in the inset data collapse has been shown
as $\frac{(1+\tanh(\kappa(x-\epsilon))}{2}$, $n=3.5$, $\kappa=0.47$ and
$\epsilon=6.31$.

It is interesting to study the delay of arrival of the signal at the receiving
end by analyzing the paths of the Active cells for the cases when
$p_{switch}=0.0$. For this we have calculated the Tortuosity $\tau$ [5], which
is defined as,

$\tau=\frac{\braket{l_{0}}}{L}=\frac{x_{0}}{L}$ (5)

where $\braket{l_{0}}=x_{0}$ is the average path length.

Of course, $\tau$ depends on percolation probability $p$, since below critical
percolation $p_{c}$ signal cannot be transmitted through the system. The
variation of $x_{0}$ as a function of $p$ is shown in Fig. 6 and $\tau$ as a
function of $p$ is shown in the inset of the same Figure in log scale for
$L=20,30,50,70$ and $100$.It has been observed that

$\tau\propto(p-p_{c})^{-u}$ (6)

The exponent $u$ is almost independent of $L$ and the approximated value is
$u\simeq 0.20$. It is to be mentioned here that value of this exponent remains
unaltered with variation of $p_{act}$.

Figure 6: $x_{0}$ vs $p$ with $p_{act}=1.00$ and $p_{switch}=0.0$ and
$L=20,30,50,70$ and $100$ n the inset data collapse has been shown as $\tau$
vs $p$ with $u\simeq 0.20$.

###  .3 Study of Cluster Distribution and Critical Exponents

Now we would move on to find out some other relevant exponents for the system.
It is well known that the cluster distribution for site percolation shows
power law behaviour for $p=p_{c}$ as $n(s)_{p_{c}}\propto s^{-{\Delta}}$ where
$\Delta$ is the Fisher exponent [11].

Figure 7: $n(s)$ vs $s$ with $L=500$, $p_{act}=1.00$ and $p_{switch}=0.0$ here
$\Delta\simeq 1.96$. Figure 8: $n(s)$ vs $s$ with $L=500$, $p_{act}=0.90$ and
$p_{switch}=0.0$ here $\Delta\simeq 1.91$.

The cluster distribution for $p_{act}=1.0$ and $p_{switch}=0.0$ for the
cluster of Inactive sites is shown in Fig. 7. Here HK algorithm has been used
to study [30] the cluster distribution. Here $\Delta$ is found to be close to
$1.9$. The distribution is shown not only for $p_{c}$ but also above and below
$p_{c}$. This nicely agrees with [31]. In Fig. 8 the cluster distribution is
shown for $p_{act}=0.9$ and $p_{switch}=0.0$.

Figure 9: $M$ vs $L$ with $p_{act}=1.00$ and $p_{switch}=0.0$ here
$d_{f}=1.81\pm 0.01$. Figure 10: Study of critical exponents, (a)
$P_{\infty}L^{\xi}$ vs $p-p_{c}$ with $\beta=0.39\pm 0.04$ and $\xi=1.95$, (b)
$\braket{s}$ vs $p-p_{c}$ with $\gamma=0.60\pm 0.02$ and (c)
$\frac{s_{max}}{L}$ vs $|p-p_{c}|$ with $\sigma=1.00\pm 0.33$, which clearly
satisfy the relation $\beta+\gamma=\frac{1}{\sigma}$.

Mass of the largest cluster $M$ for $p=p_{c}$ shows a typical behaviours
$M\sim L^{d_{f}}$ where $L$ is the system size and $d_{f}$ is the fractal
dimension. This is shown in Fig. 9. Here $d_{f}$ is found to be $d_{f}=1.81\pm
0.01$.

From the values of system dimension $d$ and fractal dimension $d_{f}$ we may
find $\Delta$ from the hyperscaling relation $\Delta=\frac{d}{d_{f}}+1$. The
value of $\Delta$ thus obtained is close to $2.10$. However, as mentioned
earlier, we have obtained a smaller value for $\Delta$.

The strength of the infinite cluster $P_{\infty}$ and the divergence of the
mean cluster size $\langle s\rangle$ are another two important quantities.
Those result in two more exponents $\beta$ and $\gamma$ which for our case
happens to be $0.39\pm 0.04$ and $0.60\pm 0.02$ and are respectively shown in
Fig. 10 (a) and (b).

From Fig. 7 another quantity is studied very often which is the cut-off
cluster size $s_{max}$. This can be found from Fig. 7 (c) by using the well
known form of $n_{s}\sim s^{-\Delta}\exp(-s/s_{max})$ for large $s$. The
corresponding variation of cut-off cluster size $s_{max}$ for our system is
shown in Fig. 10. From this we found the exponent $\sigma$ as $1.00\pm 0.33$,
the relation being $s_{max}\propto(p-p_{c})^{-1/\sigma}$.

The value of the exponent $\alpha$ can be found out easily as
$\sum_{0}^{\infty}n_{s}=(p-p_{c})^{2-\alpha}$. The value of $\alpha$ obtained
in our case is close to $2.00$.

It is clear from the above facts that the Rushbrooke’s inequality
$\alpha+2\beta+\gamma\geq 2$ is satisfied. Also the relation
$\beta+\gamma=\frac{1}{\sigma}$ is verified here to be true.

## Summary and Discussions

In summary, the work may be considered to have two aspects. Firstly, we tried
to model the heart based on percolation theory and have studied several
related features. Secondly, we have extracted the relevant exponents from the
model so obtained.

In the first part, as mentioned, the signal propagation through the heart have
been studied. It is known that the structure as well as the physiological
processes in the heart are the key factors depending on which the signal
propagates. We treat the heart as a 2D square grid of $L\times L$ cells which
may be “Waiting (0)”, “Active(1)” and “Inactive(2)”. We started from some
initial distribution of $0$s and $2$s. As the signal (AP) propagates, there
are $1$s (Active sites) too. We considered two varieties: 1. When
$2\rightarrow 0$ transition is allowed, i.e., $p_{switch}\neq 0$ and
$p_{act}=1$; 2. When $2\rightarrow 0$ transition is not allowed, i.e.,
$p_{switch}=0$ and $p_{act}\neq 1$. We studied several important quantities.
The number of arrival $N_{A}$ at the other side of the grid shows a Gaussian
behaviour when $p_{act}$ is varying and $p_{switch}=0.0$. This indicates a
kind of inhibitory behaviour and the AP/signal simply passes down the system
within a fixed time. After that no new signal can pass down through the system
as the cells are in their refractory period for an infinitely long time,
indicating a damaged heart. However, the number of arrival $N_{A}$ at the
other side of the grid shows a damped oscillatory behaviour when $p_{act}=1.0$
and $p_{switch}$ is varying, saturating to a finite value for high
$p_{switch}$. This indicates that there are always a finite number of active
cells and therefore new APs can pass through the system. It is to be noted
that very high value of $p_{switch}$ has a disadvantage as the AP is trapped
within the system which may be a cause for irregular heart rhythms.

In the other part of the work we have studied the percolation model of the
system. It has been observed that the critical probability for percolation
$p_{c}$ is close to $0.5923$ for the $p_{switch}=0$ and $p_{act}=1$. The value
of $p_{c}$ shifts towards higher and higher side as $p_{act}$ is made smaller
and smaller. Keeping $p_{act}$ fixed at $1$ and changing $p_{switch}$ a little
has severe effect on $p_{c}$ and it shifts abruptly towards $0$. The
tortuosity is studied and the related exponent is obtained as close to $0.2$.
We also focused on the study of the exponents relevant to percolation theory.
In this part we considered the parameters $p_{switch}=0$ and $p_{act}=1$ only.
The cluster distribution is studied here as shown in Fig. 7. From this we
checked the exponents $\alpha,\beta,\gamma,\sigma$ etc. The Rushbrooke
Inequality $\alpha+2\beta+\gamma\geq 2$ and the relation
$\beta+\gamma=1/\sigma$ has been observed to be satisfied.

Acknowledgement

AM acknowledges financial support from CSIR (SRF Grant no.
08/0463(12870)/2021-EMR-I). AM and SG acknowledges the computational facility
of Vidyasagar College.

## References

  * [1] S. R. Broadbent and J. M. Hammersley, Percolation processes, Mathematical Proceedings of the Cambridge Philosophical Society. 53 (3) (1957) 629–641. https://doi.org/10.1017/S0305004100032680.
  * [2] U. Dobramys, M. Mobilia, M. Pleimling and U. C. Täuber, Stochastic population dynamics in spatially extended predator–prey systems, J. Phys. A Math. Theor. 51 (2018) 063001. https://doi.org/10.1088/1751-8121/aa95c7.
  * [3] H. Barghathi, T. Vojta and J. A. Hoyos, Contact process with temporal disorder, Phys. Rev. E 94 (2016) 022111. https://doi.org/10.1103/PhysRevE.94.022111.
  * [4] E. Vigmond, A. Pashaei, S. Amraoui, H. Cochet and M. Hassaguerre, Percolation as a mechanism to explain atrial fractionated electrograms and reentry in a fibrosis model based on imaging data, Heart Rhythm 13 (2016) 1536. https://doi.org/10.1016/j.hrthm.2016.03.019.
  * [5] S. A. Niederer, J. Lumens and N. A. Trayanova, Computational models in cardiology, Nat. Rev. Cardiol. 16 (2019) 100. https://doi.org/10.1038/s41569-018-0104-y.
  * [6] R. Rabinovitch, Y. Biton, D. Braunstein, I.Aviram, R. Thieberger and A. Rabinovitch, Percolation and tortuosity in heart-like cells, Scientific Reports 11(1) (2021) 11441. https://doi.org/10.1038/s41598-021-90892-2.
  * [7] S. Biswas, S. Roy and P. Ray, Nucleation versus percolation: Scaling criterion for failure in disordered solids, Phys. Rev. E 91 (2015) 050105(R). https://doi.org/10.1103/PhysRevE.91.050105.
  * [8] A. K. Chandra Percolation in a kinetic opinion exchange model, Physical Review E 85 (2012) 021149. https://doi.org/10.1103/PhysRevE.85.021149.
  * [9] F. Bagnoli, G. de B. Cavalcabo’, A Simple Model of Knowledge Percolation, archive, arXiv:2206.03267 (2022). https://doi.org/10.48550/arXiv.2206.03267.; B. Chopard, S. Bandini, A. Dennunzio, A. Haddad, Mira (eds) Cellular Automata. ACRI 2022. Lecture Notes in Computer Science, Springer, Cham. 13402 (2022). https://doi.org/10.1007/978-3-031-14926-9.
  * [10] P. Bak, K. Chen and C. Tang, A forest-fire model and some thoughts on turbulence, Phys. Lett. A 147 (1990) 297. http://dx.doi.org/10.1016/0375-9601(90)90451-S.; P. Bak, C. Tang and K. Wiesenfeld, Self-organized criticality, Phys. Rev. A 38 (1988) 364. https://doi.org/10.1103/PhysRevA.38.364.
  * [11] D. Stauffer and A. Aharony, Introduction to Percolation Theory, Taylor & Francis, (1994). https://doi.org/10.1201/9781315274386.
  * [12] H. S. Song and S. H. Lee, Effects of wind and tree density on forest fire patterns in a mixed-tree species forest, Forest Science and Technology, 13(1) (2016) 9-16. http://dx.doi.org/10.1080/21580103.2016.1262793.
  * [13] B. Drossel, and F. Schwabl, Self-organized critical forest-fire model, Phys. Rev. Lett. 69 (1992) 1629. https://doi.org/10.1103/PhysRevLett.69.1629.
  * [14] W. K. Mo${\beta}$ner,B. Drossel, and F. Schwabl, Computer simulations of the forest-fire model, Phys. A Stat. Mech. Appl. 190(3) (1992) 205-217. https://doi.org/10.1016/0378-4371(92)90032-L.
  * [15] D. Knežević, M. Knežević. Semi-directed percolation in two dimensions, Physica A 444 (2016) 560–565. https://doi.org/10.1016/j.physa.2015.10.079.
  * [16] S. Redner, Conductivity of random resistor-diode networks, Phys. Rev. B 25 (1982) 5646–5655. https://doi.org/10.1103/PhysRevB.25.5646.
  * [17] H. O. Mártin and J. Vannimenus, Partially directed site percolation on the square and triangular lattices, J. Phys. A: Math. Gen. 18 (1985) 1475-1482. 10.1088/0305-4470/18/9/027.
  * [18] R. F. Gibson, A review of recent research on mechanics of multifunctional composite materials and structures, Composite Structures 92, (2010) 2793–2810. http://dx.doi.org/10.1016/j.compstruct.2010.05.003.
  * [19] V. Langlois, V. H. Trinh, and C. Perrot, Electrical conductivity and tortuosity of solid foam: effect of pore connections, Phys. Rev. E 100 (2019) 013115. https://doi.org/10.1103/PhysRevE.100.013115.
  * [20] B. Ghanbarian, A. G. Hunt, M. Sahimi, R. P. Ewing, and T. E. Skinner, Percolation theory generates a physically based description of tortuosity in saturated and unsaturated porous media, Soil Sci. Soc. Am. J. 77 (2013) 1920. http://dx.doi.org/10.2136/sssaj2013.01.0089.
  * [21] D. S. Desai and S. Hajouli, Arrhythmias, StatPearls [Internet]. Treasure Island (FL): StatPearls Publishing, (2023) PMID: 32644349.
  * [22] O. Dössel , M. W. Krueger, F. M. Weber, M. Wilhelm and G. Seemann, Computational modeling of the human atrial anatomy and electrophysiology, Med. Biol. Eng. Comput. 50 (2012) 773–799. https://doi.org/10.1007/s11517-012-0924-6.
  * [23] A. Shirakabe, Y. Ikeda, S. Sciarretta, D. K. Zablocki, J. Sadoshima, Aging and Autophagy in the Heart. Circ Res. 118(10) (2016) 1563-1576, PMID: 27174950, PMCID: PMC4869999. https://doi.org/10.1161%2FCIRCRESAHA.116.307474.
  * [24] P. Anversa, A. Leri, Innate regeneration in the aging heart: healing from within, Mayo Clin Proc. 88(8) (2013) 871-883, PMID: 23910414, PMCID: PMC3936323. https://doi.org/10.1016/j.mayocp.2013.04.001
  * [25] D. L. Johnson, J. Koplik and R. Dashen, Theory of dynamic permeability and tortuosity in fluidsaturated porous media, Journal of Fluid Mechanics, 176 (1987) 379402. https://doi.org/10.1017/S0022112087000727.
  * [26] B. Derrida and D. Stauffer, Corrections to scaling and phenomenological renormalization for 2-dimensional percolation and lattice animal problems, Journal de Physique 46 (1985) 1623. http://dx.doi.org/10.1051/jphys:0198500460100162300.
  * [27] K. Malarz and S. Galam, Square-lattice site percolation at increasing ranges of neighbor bonds, Physical Review E 71 (2005) 016125. https://doi.org/10.1103/PhysRevE.71.016125.
  * [28] X. Feng, Y. Deng, H. W. J. Blöte, Percolation transitions in two dimensions, Physical Review E 78 (2008) 031136. https://doi.org/10.1103/PhysRevE.78.031136.
  * [29] S. Mertens, Exact site-percolation probability on the square lattice, J. Phys. A: Math. Theor. 55 (2022) 334002. https://doi.org/10.1088/1751-8121/ac4195.
  * [30] J. Hoshen and R. Kopelman, Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm, Phys. Rev. B 14 (1975) 3438. https://doi.org/10.1103/PhysRevB.14.3438.
  * [31] M. Nicodemi, Percolation, Complexity Science & Theor. Phys., University of Warwick. lesson percolation.pdf


# Correcting for Reporting Delays in Cyber Incidents

Seema Sangari School of Data Science and Analytics, Kennesaw State University,
3391 Town Point Dr. NW,Kennesaw, GA 30144. Contact:
<EMAIL_ADDRESS>Dr. Eric Dallal AIR Worldwide, Verisk Cyber
Solutions, Lafayette City Center, 2 Avenue de Lafayette, 2nd Floor, Boston, MA
02111. Contact<EMAIL_ADDRESS>
###### Abstract

With an ever evolving cyber domain, delays in reporting incidents are a well-
known problem in the cyber insurance industry. Addressing this problem is a
requisite to obtaining the true picture of cyber incident rates and to model
it appropriately. The proposed algorithm addresses this problem by creating a
model of the distribution of reporting delays and using the model to correct
reported incident counts to account for the expected proportion of incidents
that have occurred but not yet been reported. In particular, this correction
shows an increase in the number of cyber events in recent months rather than
the decline suggested by reported counts.

###### keywords:

Delay Distribution, Optimization, Reporting Delays

##  1 Introduction

The cyber security domain is evolving rapidly, with new attack vectors
emerging regularly, and even the most up-to-date data cannot be considered
complete. Cyber incidents take a long time to become known and even longer to
appear in online databases. While some cyber events are known immediately
after they occur, most events are often reported many months or years after
the event actually occurred, resulting in biased data. Marriott’s major cyber
incident occurred in 2014 but was only reported in 2018. Major cyber events
become headlines in leading newspapers when reported publicly rather than when
they happened. Smaller cyber events may never be reported at all, or have
extreme delays, as only public companies and those with personally
identifiable information may be obligated to report. Sometimes reporting can
take 5-10 years, for various intentional or unintentional reasons - failing to
realize that a cyber incident happened, failing to immediately determine the
extent of accessed or stolen data, or deciding not to publicize the incident
for fear of reputation risk and consequent financial impacts. As a result,
reporting delays are often observed in historical cyber event databases. These
databases show a decrease in cyber incidents in the recent few years. Coleman
et al. (2021) raised the concern that cyber incidents would remain undetected
due to advanced threat techniques.

Cyber risk modeling firms rely upon historical data to build their models,
which are in turn relied upon by cyber insurers for underwriting, portfolio
management, and risk transfer. To build robust loss estimation models for
today’s evolving cyber world, the most recent and updated information is
required, with as little bias as possible. Correcting reporting delays in
these databases is therefore a key requirement to have trustworthy cyber
insurance models. With the necessary corrections, one can then properly
examine temporal trends in the targeting of industries or in attacker tactics.

Terminology followed in this paper:
$i$: Incident
$I$: Set of all incidents
$delay,\delta_{i}$: Time between the incident date and the reporting date
$age,A_{i}$: Time between the incident date and the last incident reporting
date in the data
$f_{\Delta}$: Probability density function of the $delay$ distribution
$F_{\Delta}$: Cumulative distribution function of the $delay$ distribution

##  2 Literature Review

Most of the literature on reporting delays addresses the medical space. We are
not aware of any existing papers that propose a methodology for correcting
cyber incident counts to account for reporting delays. However, Coleman et al.
(2021) does examine both the distribution of the number of days to discover
cyber incidents and the number of days to disclose them.

Harris (1987) described reporting delays as a statistical problem for the
first time. Heisterkamp et al. (1988a, b, 1989) and Brookmeyer and Damiano
(1989) made distributional assumptions and built linear/quadratic models
whereas Rosenberg (1990) and Cheng and Ford (1991) suggested Poisson models.
These are easy to implement but assume stationary reporting delays and don’t
capture trends. Morgan and Curran (1986), Downs et al. (1987, 1988)), Healy
and Tillett (1988) and Heisterkamp et al. (1988a, b, 1989) fitted exponential,
integrated logistic and log-linear models to capture trends. Gail and
Brookmeyer (1988), Brookmeyer and Liao (1990), Kalbfleisch and Lawless (1991)
and Esbjerg et al. (1999) applied conditional probabilities to capture trends
but this resulted in over-fitting. Lawless (1994) proposed a multinomial model
with distributed random effects based on Dirichlet/Poisson/Gamma distributions
to capture trends in a timely fashion but failed to handle longer delays.

Wang (1992) suggested maximum likelihood estimation (MLE) based non-parametric
and semi-parametric approaches but with complete111Complete data - No further
events are expected to be reported with delays. data. Harris (2020) suggested
correcting COVID cases with an expectation-maximization (EM) algorithm and
trained the model with complete data to correct test data. White et al. (2009)
and Weinberger et al. (2020) proposed a simpler method based on proportions
but also require complete data to train.

Wang (1992), Keiding and Moeschberger (1992) and Midthune et al. (2005)
applied truncated models to avoid random effects but require stable reporting
delays.

Höhle and An Der Heiden (2014), Bastos et al. (2019) and Chitwood et al.
(2020) suggested a Bayesian and hierarchical approach with Poisson and
negative binomial distributions. It is easy to implement but makes
distributional assumptions. Noufaily et al. (2015, 2016)) suggested a log-
likelihood approach with a truncation model that is data driven but sensitive
to the choice of three reporting time steps.

Bastos et al. (2019) suggested a chain-ladder approach but it is sensitive to
outliers.

Jewell (1989), Zhao et al. (2009), Zhao and Zhou (2010), and Avanzi et al.
(2016) investigated cyber claims data to account for reporting delays from a
capital reserving perspective. This problem is different from the one being
investigated, since reporting delays in claims are due only to detection
delays.

As stated in Brookmeyer and Liao (1990), none of these approaches deal with
delays longer than any previously observed.

##  3 Data

The data used is a proprietary data set constructed by merging multiple source
cyber event sets together. The source cyber event sets contained differing
fields that, at minimum, provided information as to the “what”, “when” and
“who” of an event. Concretely, the data sets included a description of the
event, the name of the company to which the cyber event occurred (N.B.:
aggregation events separately list each company known to be impacted), along
with the date that the event occurred and the date that the event was
reported.

Because company names frequently differ from one data set to another, the data
sets could not be de-duplicated based on direct string matching. Instead, the
event data sets were matched to a firmographic data set containing
approximately 55 million businesses in the US. This was done via a previously
developed matching algorithm that examined company name, industry
classification (e.g., via NAICS codes), address information, and any other
fields common to both the cyber event data set and the firmographic data set.
Events in distinct cyber event data sets were identified as the same when:

  1. 1.

They were matched to the same company in the firmographic data set; and

  2. 2.

They were listed as having occurred within 1 week of each other.

Limitations: Some events do not have an occurrence date listed, and were
therefore excluded from this analysis. There are also large spikes in event
counts listed as having occurred on January 1st in most years (Fig. 1), which
was assumed to be a default value when only the year of the event was known.
These events were therefore re-distributed proportionally throughout the year
as shown in Fig. 1 but excluded while developing the approach.

With Default Date

Adjusted for Default Date

Figure 1: Cyber Event Counts until December 2018

##  4 Proposed Approach

The proposed approach consists of estimating the reporting delay distribution
$f_{\Delta}$ from the empirical data. With this distribution, a corrected
count of events with age $a$ can be determined by dividing the raw counts by
$F_{\Delta}(a)$, as the latter expression represents the proportion of events
that are reported within a delay, $\delta$, of less than $a$ or, equivalently,
the proportion of events that are reported as of the present. There are four
complications with estimating the delay distribution in the “obvious” way:

  * •

The nature of reporting delays means that direct estimates from empirical data
will be biased towards shorter delays, since recent events could only appear
in the data set in the first place if the reporting delay is small.

  * •

The estimates assume zero probability of any delay longer than the longest in
the data set.

  * •

The estimates for long reporting delays are based on few data points.

  * •

The reporting delay distribution may not be stationary.

Mathematically, the first problem is that, for incidents with age $a$, the
delay $\delta_{i}$ can only be observed if $\delta_{i}<=a$. This problem can
be resolved if one assumes that $f_{\Delta}$ is stationary, in which case the
distribution can be estimated for delays, $\delta$, larger than $a$ from older
events. This is the approach taken in Algorithm 1. To test for stationary, the
algorithm was run on two year windows of the data set centered at different
times. (Note that the two year windows include all events that _occur_ within
the window, irrespective of when they are reported.) This showed that the
delay distribution is non-stationary, as can be seen in Fig. 2.

Probability Distribution

Cumulative Distribution

Figure 2: PDF and CDF of Delay Distribution generated for Dec.,’12 and
Dec.,’16

To deal with this non-stationarity, the delay distribution was modeled via
parametric distributions where parameters were estimated over monthly rolling
two year windows. For each two year window, optimal parameters were determined
based on matching cumulative distributions. Specifically, Algorithm 1 was run
on each of the monthly two year windows to obtain an estimate of the delay
distribution for that time window. For a given two year window and modeled
distribution, an optimization algorithm was used to determine the parameters
that gave the best match between:

  * •

the delay distribution obtained by Algorithm 1; and

  * •

the parametric distribution restricted to the domain where the maximum delay
is less than the window’s maximum reporting delay.

The optimal distribution parameters for each window were computed using a
derivative free optimizer.

In what follows, the above two distributions will be referred to as the
“debiased empirical delay distribution” and the “modeled delay distribution”,
respectively.

###  4.1 Generating the Debiased Empirical Delay Distribution

Inspired from Brookmeyer and Liao (1990), the proposed algorithm works on the
limitation that the delay, $\delta$, is not considered beyond age and hence
only a distribution conditioned on the delay being less than or equal to age
can be estimated. The algorithm corrects the incident counts with $F_{\Delta}$
estimated from empirical data. The algorithm applies a top down approach to
estimate the distribution “from the outside in”, accounting for the estimated
proportion of events that have occurred but not yet been reported as the
distribution itself is being computed.

Let $A_{\max}:=\max_{i\in\mathcal{I}}A_{i}$ be the maximal age of any event in
the data set. Also let $h_{A}(a)$ be the number of incidents of age $a$, and
let $h_{\Delta}(\delta)$ be the number of incidents with delay $\delta$.
Formally,

$\displaystyle h_{A}(a)$ $\displaystyle=$
$\displaystyle|\\{i\in\mathcal{I}:A_{i}=a\\}|$ (1) $\displaystyle
h_{\Delta}(\delta)$ $\displaystyle=$
$\displaystyle|\\{i\in\mathcal{I}:\Delta_{i}=\delta\\}|$ (2)

Then the equation to estimate the delay distribution is:

$\displaystyle
f_{\Delta}(\delta)=\frac{h_{\Delta}(\delta)}{\sum_{a=\delta}^{A_{\max}}\nicefrac{{h_{A}(a)}}{{F_{\Delta}(a)}}}$
(3)

Intuitively, the distribution is generated based on the ratio of the number of
events with the given delay period, $h_{\Delta}(\delta)$, to the best estimate
of the true number of events whose age is old enough to be seen in the
incident data set (i.e., where the age is larger than the delay under
consideration),
$\sum_{a=\delta}^{A_{\max}}\nicefrac{{h_{A}(a)}}{{F_{\Delta}(a)}}$. Algorithm
1 shows the algorithmic implementation to generate the distribution.

1:The histograms, $h_{A}$ and $h_{\Delta}$, computed as in Eqs. (1) and (2),
respectively.

2:The distribution $f_{\Delta}$.

3:

4:function ComputeDelayDistribution($h_{A}$, $h_{\Delta}$)

5: $A_{\max}\leftarrow\max_{i\in\mathcal{I}}A_{i}$

6: $F_{\Delta}(A_{\max})\leftarrow 1$

7: $\delta_{max}\leftarrow\max_{i\in\mathcal{I}}\delta_{i}$

8: for $\delta\leftarrow A_{\max}$ to $\delta=0$ do

9: $den\leftarrow 0$

10: for $a\leftarrow\delta$ to $\delta_{\max}$ do

11: $den\leftarrow den+h_{A}(a)/F_{\Delta}(a)$ $\triangleright$ Computes
denominator

12: end for

13: $f_{\Delta}(\delta)\leftarrow h_{\Delta}(\delta)/den$ $\triangleright$
Computes PDF

14: $F_{\Delta}(\delta-1)\leftarrow F_{\Delta}(\delta)-f_{\Delta}(\delta)$
$\triangleright$ Updates CDF

15: $\delta_{max}\leftarrow\delta$

16: end for

17: return $f_{\Delta}$

18:end function

Algorithm 1 Algorithm for computing the debiased empirical delay distribution

###  4.2 Generating the Modeled Delay Distribution

As previously described, the modeled delay distribution is determined by an
optimization that minimizes the difference in the cumulative distribution
functions of the debiased empirical delay distribution and the modeled delay
distribution, restricted to the domain $[0,\delta_{\max}]$ (N.B.:
$\delta_{\max}$ is the maximum observed delay). The PDF and CDF plots shown in
Fig. 2 suggest that a single distribution will not provide a good fit, as the
delay distribution is bi-modal. Rather, a mixture of distributions would be
required. Fig. 2 suggests a mix of exponential and normal distributions.

Mathematically,

$\displaystyle F_{\theta}(\delta)$ $\displaystyle=\alpha
F_{Exp}(\delta,Scale)+(1-\alpha)F_{N}(\delta,\mu,\sigma)$ (4)
$\displaystyle\text{where }F_{Exp}=\text{Exponential CDF with parameter,
$Scale=\nicefrac{{1}}{{\lambda}}$}$ $\displaystyle\quad\quad\quad
F_{N}=\text{Normal CDF with parameters, $\mu$ and $\sigma$}$

A possible interpretation of the bi-modal nature of the reporting delay
distribution is that there are two distributions: one for events that are
discovered almost immediately (modeled by the exponential) and one for events
which have a delay due to both discovery time and public disclosure time
(modeled by the normal). The parameter $\alpha$ can therefore be interpreted
as the proportion of events that are discovered right away by the
organization.

Since the normal distribution is defined on $(-\infty,\infty)$ and reporting
delays cannot be negative, the modeled delay distribution needs to be
adjusted. The CDF in Eq. 4 truncated to $[0,\infty)$ (and re-normalized) can
be expressed as

$\displaystyle F_{\theta}(\delta)=$
$\displaystyle\frac{\alpha(F_{Exp}(\delta,Scale))+(1-\alpha)\overbrace{(F_{N}(\delta,\mu,\sigma)-F_{N}(0,\mu,\sigma)}^{\text{Truncated
Normal Distribution until
$\delta$}}}{\alpha+(1-\alpha)\underbrace{(1-F_{N}(0,\mu,\sigma)}_{\text{Truncated
Normal Distribution over $[0,\infty)$}}}$ (5) where $0\leq\delta\leq\infty$

Since the debiased empirical delay distribution is only defined on
$[0,\delta_{\max}]$, which is the domain on which it is compared to the
modeled delay distribution, the truncation of the modeled distribution to this
domain is also defined, denoted by $F^{\prime}_{\theta}$ as below.

$\displaystyle F^{\prime}_{\theta}(\delta)=$
$\displaystyle\frac{\alpha(F_{Exp}(\delta,Scale))+(1-\alpha)\overbrace{(F_{N}(\delta,\mu,\sigma)-F_{N}(0,\mu,\sigma)}^{\text{Truncated
Normal Distribution until
$\delta$}}}{\alpha(F_{Exp}(\delta_{\max},Scale))+(1-\alpha)\underbrace{(F_{N}(\delta_{\max},\mu,\sigma)-F_{N}(0,\mu,\sigma)}_{\text{Truncated
Normal Distribution over $[0,\delta_{\max}]$}}}$ (6) where
$0\leq\delta\leq\delta_{\max}$

###  4.3 Defining the Optimization Function

Defining the optimization function was challenging due to two factors in
particular:

  * •

There are many combinations of parameters that give approximately the same
distribution when restricted to the domain $[0,\delta_{\max}]$, but which
differ substantially in how much of the total distribution’s weight is
contained in this domain.

  * •

As two year windows closer to the present are considered, the quantity of data
shrinks, resulting in increasingly unstable parameter estimates.

The optimization function used is shown in Eq. 7 below. The first
term222$F^{\prime}_{\theta}$ computed as defined in Eq. 6
$\|\log_{10}F^{\prime}_{\theta}-\log_{10}F_{\Delta}\|^{2}$ reduces the CDF
difference between the debiased empirical delay distribution and the modeled
delay distribution over the domain $[0,\delta_{\max}]$. The purpose of
applying $\log_{10}$ weights is to place more emphasis on a good fit for the
initial months333 The rationale behind $\log_{10}$, is to obtain CDF values
close from the point of the ratio between the two distributions (modeled and
debiased empirical), not in terms of absolute difference - a $\log_{10}$ CDF
difference between 0.03 and 0.06 would be a factor of two difference in the
correction whereas a difference between 0.93 and 0.96 would be much smaller
despite the fact that the absolute difference is 0.03 in both cases.. As
mentioned above, there are many different combinations of parameters that
would result in comparable errors in the first optimization term, but which
nevertheless differ substantially over the domain $[0,\infty)$. This problem
is relatively minor when the range $[0,\delta_{\max}]$ contains the bulk of
the distribution, which it does when $\delta_{\max}$ is substantially greater
than the second peak of the debiased empirical delay distribution. But as two
year windows closer to the present are considered, this ceases to be the case.
In order to avoid this problem, the optimization function must consider
modeled delay distribution values beyond $\delta_{\max}$. The second term
$\|\log_{10}S_{\theta}-\log_{10}S_{\theta^{\prime}}\|^{2}$ penalizes large
differences between the CDFs of consecutive modeled distributions beyond
$\delta_{\max}$. In order to further minimize parameter instability for recent
two year windows, another set of weights is assigned to the first two terms,
which effectively diminishes the importance of a good fit with the empirical
data and increases the importance of parameter stability as the amount of
available data diminishes.

The third and fourth terms are penalization terms - $F_{N}^{2}(0,\mu,\sigma)$
term penalizes negative delays introduced by the normal distribution (defined
over $(-\infty,+\infty)$) whereas $S_{\theta}^{2}(10Y)$ term penalizes delays
beyond 10 years.

Mathematically, the optimization function is defined as

$\displaystyle\theta_{Opt}$
$\displaystyle=\underset{\theta=(\alpha,Scale,\mu,\sigma)}{argmin}\frac{\delta_{\max}}{\delta_{Fix}}\underbrace{\|\log_{10}F^{\prime}_{\theta}-\log_{10}F_{\Delta}\|^{2}}_{\text{$\delta\in[0,\delta_{\max}]$}}$
(7)
$\displaystyle\qquad\qquad\qquad+\left(1-\frac{\delta_{\max}}{\delta_{Fix}}\right)\underbrace{\|\log_{10}S_{\theta^{\prime}}-\log_{10}S_{\theta}\|^{2}}_{\text{$\delta\in(\delta_{\max},\delta_{Fix}]$}}$
$\displaystyle\qquad\qquad\qquad+\underbrace{F_{N}^{2}(0,\mu,\sigma)}_{\text{$\delta<0$}}+\underbrace{S_{\theta}^{2}(10Y)}_{\text{$\delta>10Years$}}$
where $\delta_{Fix}$ is the Maximum value of $\delta$ in the dataset.

$F^{\prime}_{\theta}$ is as defined in Eq. 6 whereas $S_{\theta}$ is the
survival function, defined as the complement of $F_{\theta}$ of Eq. 5:

$\displaystyle S_{\theta}$ $\displaystyle=1-F_{\theta}$ (8)

Finally, $\theta^{\prime}$ refers to the previous two year window’s optimal
parameters. Since no previous parameters are available for the first two year
window, the second term is taken to be zero for that window.

$\displaystyle\|\log_{10}S_{\theta^{\prime}}-\log_{10}S_{\theta}\|^{2}=0$ (9)
$\theta^{\prime}$ refers to optimal parameters at previous step.

The covariance matrix adaptation evolution strategy (CMA-ES) was applied to
compute the modeled distribution parameters. It is a derivative free
optimization algorithm, a type typically used when derivatives are difficult
or costly to compute (Hansen (2006, 2016, 2019)).

###  4.4 Computing the Corrected Counts

Once the modeled delay distributions (one for each two year window) have been
obtained by optimization, event counts can be corrected based on the
cumulative distribution function of the full modeled distribution (i.e.,
defined over $[0,\infty)$) defined by Eq. 5.

Specifically,

Corrected Count $\displaystyle=\frac{\text{Reported Counts for month,
`m'}}{F_{\theta}(a)}$ (10) where $a$ is age of the given month, ‘m’.

From each modeled distribution, the correction factor is computed at only one
point, $age$, of the modeled distribution for the given month’s correction.

##  5 Results

Fig. 3 shows two examples of plots comparing PDFs of the fitted parametric
modeled distribution, its truncation to the domain $[0,\delta_{\max}]$ — where
$\delta_{max}$ is computed for each window individually, and the debiased
empirical delay distribution. Fig. 3 shows this comparison for the two year
window starting from July 2012 until June 2016 and Fig. 3 shows this
comparison for the most recent window starting from January, 2017 until
December 2018.

###  5.1 Parameters and Interpretation

Fig. 4 shows the parameter plots of the delay distribution generated for each
monthly two year rolling window.

From July 2012 to June, 2014

From January, 2017 to December, 2018

Figure 3: Comparing PDFs of Debiased Delay Distribution with Parametric
Modeled Distribution

The alpha plot (Fig. 4) suggests that organizations discover 8-18% of cyber
events right away ($8\%\leq\alpha\leq 18\%$).

The scale plot (Fig. 4) suggests that the short delays modeled by the
exponential distribution had a mean of less than 60 days delay until early
2016 but increased rapidly to around 140 days in early 2018.

The normal distribution mean, $\mu$, (Fig. 4) and standard deviation,
$\sigma$, (Fig. 4) parameter plots suggest that the longer delays modeled by
the normal distribution remained consistent over time. The period of longer
delays remain consistent varying within $\pm 10\%$ range.

Alpha Parameter

Scale Parameter

Mean Parameter

Sigma Parameter

Figure 4: Plots of Modeled Distribution Parameters based on Empirical Debiased
Delay Distribution

###  5.2 Corrections and Validation

Fig. 5 shows the corrected incident counts based on the proposed methodology.
Figs. 5 and 5 show the corrected counts for the events reported by Dec. 2017
and by Dec. 2018, respectively. Although the corrections follow similar trends
in both, the correction factors vary substantially.

To validate the proposed algorithm, the counts reported until December 2017
(2018) were corrected for a year ahead and compared against the counts
reported as of December 2018 (2019). The year ahead correction factor is
computed as

$\text{ Year ahead
}F_{\theta}(a,a+1\>Year)=\frac{F_{\theta}(a)}{F_{\theta}(a+1\>Year)}$ (11)

where $a$ is the age of the event counts being corrected.

Whereas the 2017 year ahead corrections (Fig. 5) initially show close
agreement with the 2018 counts, more recent year ahead corrections
underestimate the 2018 counts. On the other hand, the 2018 year ahead
corrections (Fig. 5) generally overestimate the 2019 counts, except for the
most recent months, which show close agreement. As stated in section 4.3, the
debiased empirical delay distribution has fewer data points for more recent
two year windows so weights of ($\nicefrac{{\delta_{\max}}}{{\delta_{Fix}}})$
and ($1-\nicefrac{{\delta_{\max}}}{{\delta_{Fix}}}$) are used in the
optimization function to dynamically adjust the weight given to the CDF before
and after $\delta_{\max}$ respectively. By removing these weights, better
estimates for recent months might be obtained but would come at the cost of
more parameter instability and worse validation plots (overfitting). In either
Fig. 5 or Fig. 5, the corrected counts (dashed line) show a trend of
increasing incident counts since 2016, which is contrary to the diminishing
trend seen in the raw counts. The trend in corrected counts is therefore much
more in line with reports from insurers and other organizations that release
reports on cyber risk.

2017 corrections Vs 2018 cumulative counts

2018 corrections Vs 2019 cumulative counts

Until 201X - Counts reported as of 201X adjusted for the default date of
January 1, proportionally 201X Corrected - “Counts until 201X” corrected based
on Eq. 10 201X Corrected Year Ahead - “Counts until 201X” corrected based on
Eq. 11

Figure 5: Validation Plots

##  6 Conclusion

This work examined the long known problem of reporting delays in historical
cyber events databases and proposed an algorithm to correct for these delays.
Interestingly, the true distribution of reporting delays appears to be bi-
modal, which we have interpreted as a mixture of two distributions: one for
incidents that are discovered immediately, modeled by an exponential
distribution, and one for incidents that are not immediately discovered,
modeled by a normal distribution. With this form of reporting delay
distribution, we obtained non-stationary modeled delay distributions via
optimization. These modeled delay distributions were used to estimate the
total number of cyber incidents that will eventually be reported from the
current counts. The approach was validated by estimating year ahead
corrections.

To understand the current cyber threat landscape and to create robust cyber
risk models, one needs accurate historical data. While it is not possible to
get the exact count of cyber events, the proposed algorithm aims to correct
for reporting delays approximately. The reported cyber incident counts in
recent times show a decreasing trend simply because incidents have not been
reported yet, even though they have actually already occurred. However, in
reality, the rate of cyber incidents is increasing and that is what the
algorithm reveals.

The general approach can be applied to any form of data with reporting delays
including other long tailed insurance claims (such as liability), and COVID-19
pandemic data. The current study corrects the overall counts for US cyber
events, and the industry specific correction of cyber event counts is a future
direction of this work.

## Acknowledgment

The study was done in collaboration with AIR Worldwide using their proprietary
cyber data. The authors would like to thank Scott Stransky for his comments
and suggestions that helped improve the approach.

## References

  * Avanzi et al. (2016) Benjamin Avanzi, Bernard Wong, and Xinda Yang.  A micro-level claim count model with overdispersion and reporting delays.  _Insurance: Mathematics and Economics_ , 71:1–14, 2016\.  ISSN 01676687.  10.1016/j.insmatheco.2016.07.002.
  * Bastos et al. (2019) Leonardo S. Bastos, Theodoros Economou, Marcelo F.C. Gomes, Daniel A.M. Villela, Flavio C. Coelho, Oswaldo G. Cruz, Oliver Stoner, Trevor Bailey, and Claudia T. Codeço.  A modelling approach for correcting reporting delays in disease surveillance data.  _Statistics in Medicine_ , 38(22):4363–4377, 2019\.  ISSN 10970258.  10.1002/sim.8303.
  * Brookmeyer and Damiano (1989) Ron Brookmeyer and Anne Damiano.  Statistical methods for short‐term projections of AIDS incidence.  _Statistics in Medicine_ , 8(1):23–34, 1989.  ISSN 10970258.  10.1002/sim.4780080105.
  * Brookmeyer and Gail (1986) Ron Brookmeyer and Mitchell H. Gail.  Minimum Size of the Acquired Immunodeficiency Syndrome (Aids) Epidemic in the United States.  _The Lancet_ , 328(8519):1320–1322, 1986.  ISSN 01406736.  10.1016/S0140-6736(86)91444-3.
  * Brookmeyer and Liao (1990) Ron Brookmeyer and Jiangang Liao.  The analysis of delays in disease reporting: Methods and results for the acquired immunodeficiency syndrome.  _American Journal of Epidemiology_ , 132(2):355–365, 1990.  ISSN 00029262.  10.1093/oxfordjournals.aje.a115665.
  * Cheng and Ford (1991) Fenny F.K. Cheng and Wesley L. Ford.  Adjustment of aids surveillance data for reporting delay to the editor:, 1991.  ISSN 10779450.
  * Chitwood et al. (2020) Melanie H Chitwood, Marcus Russi, Kenneth Gunasekera, Joshua Havumaki, Virginia E Pitzer, Joshua L Warren, Daniel Weinberger, Ted Cohen, and Nicolas A Menzies.  Bayesian nowcasting with adjustment for delayed and incomplete reporting to estimate COVID-19 infections in the United States.  _medRxiv_ , page 2020.06.17.20133983, 2020.  10.1101/2020.06.17.20133983.  URL http://medrxiv.org/content/early/2020/06/20/2020.06.17.20133983.abstract.
  * Coleman et al. (2021) Derryck Coleman, Madeleine Conley, and Nicole Hallas.  Trends in Cybersecurity Breaches.  Technical report, Audit Analytics, Massachusetts, USA, 2021.  URL https://go.auditanalytics.com/cybersecurityreport.
  * Downs et al. (1987) A. M. Downs, R. A. Ancelle, H. J.C. Jager, and J. B. Brunet.  AIDS in Europe: Current trends and short-term predictions estimated from surveillance data, January 1981-June 1986.  _Aids_ , 1(1):53–57, 1987.  ISSN 02699370.
  * Downs et al. (1988) A. M. Downs, R.A. Ancelle, J. C. Jager, S. H. Heisterkamp, J. A.M. Van Druten, E. J. Ruitenberg, and Brunet J.B.  The statistical estimation, from routine surveillance data, of past, present and future trends in AIDS incidence in Europe.  In J. C. Jager and E. J. Ruitenberg, editors, _Statistical Analysis and Mathematical Modelling of AIDS_ , pages 1–16. Oxford University Press, Oxford, 1988.
  * Esbjerg et al. (1999) Sille Esbjerg, Niels Keiding, and Nils Koch-Henriksen.  Reporting delay and corrected incidence of multiple sclerosis.  _Statistics in Medicine_ , 1999.  ISSN 02776715.  10.1002/(SICI)1097-0258(19990715)18:13¡1691::AID-SIM160¿3.0.CO;2-D.
  * Gail and Brookmeyer (1988) Mitchell H. Gail and Ron Brookmeyer.  Methods for projecting course of acquired immunodeficiency syndrome epidemic, 1988.  ISSN 00278874.
  * Hansen (2006) N Hansen.  The CMA evolution strategy: a comparing review, in: J.A. Lozano, P. Larranaga, I. Inza, E. Bengoetxea (Eds.), Towards A New Evolutionary Computation. Advances on Estimation of Distribution Algorithms, Springer.  192:75–102, 2006.
  * Hansen (2016) Nikolaus Hansen.  The CMA Evolution Strategy: A Tutorial.  _Computing Research Repository_ , apr 2016.  URL http://arxiv.org/abs/1604.00772.
  * Hansen (2019) Nikolaus Hansen.  CMA - Python Package, 2019.  URL https://pypi.org/project/cma/.
  * Harris (1987) Jeffrey E Harris.  Delay in Reporting Acquired Immune Deficiency Syndrome (AIDS).  _National Bureau of Economic Research Working Paper Series_ , No. 2278, 1987.  URL https://www.nber.org/papers/w2278.
  * Harris (2020) Jeffrey E Harris.  Overcoming Reporting Delays Is Critical to Timely Epidemic Monitoring: The Case of COVID-19 in New York City.  _medRxiv_ , page 2020.08.02.20159418, 2020.  10.1101/2020.08.02.20159418.  URL https://doi.org/10.1101/2020.08.02.20159418.
  * Healy and Tillett (1988) M. J. R. Healy and H. E. Tillett.  Short-Term Extrapolation of the AIDS Epidemic.  _Journal of the Royal Statistical Society. Series A (Statistics in Society)_ , 151(1):50, 1988.  ISSN 09641998.  10.2307/2982184.
  * Heisterkamp et al. (1988a) S. H. Heisterkamp, J. C. Jager, A. M. Downs, and J. A.M. Van Druten.  The use of Genstat in the estimation of expected numbers of AIDS cases adjusted for reporting delays.  In _Fifth Genstat Conference_ , pages 4–18, 1988a.
  * Heisterkamp et al. (1988b) S. H. Heisterkamp, J. C. Jager, A. M. Downs, J. A.M. Van Druten, and E. J. Ruitenberg.  Statistical estimation of AIDS incidence from surveillance data and the link with modelling of trends.  In _Statistical Analysis and Mathematical Modelling of AIDS_ , pages 17–25. Oxford University Press, Oxford, 1988b.
  * Heisterkamp et al. (1989) S. H. Heisterkamp, J. C. Jager, E. J. Ruitenberg, J. A.M. Van Druten, and A. M. Downs.  Correcting reported aids incidence: A statistical approach.  _Statistics in Medicine_ , 8(8):963–976, 1989\.  ISSN 10970258.  10.1002/sim.4780080807.
  * Höhle and An Der Heiden (2014) Michael Höhle and Matthias An Der Heiden.  Bayesian nowcasting during the STEC O104: H4 outbreak in Germany, 2011\.  _Biometrics_ , 70(4):993–1002, 2014.  ISSN 15410420.  10.1111/biom.12194.
  * Jewell (1989) WS S. Jewell.  Predicting ibnyr events and delays.  _ASTIN Bulletin_ , 19(I):25–56, 1989.  ISSN 17831350.
  * Kalbfleisch and Lawless (1991) J D Kalbfleisch and J F Lawless.  Regression Models for Right Truncated Data With Applications To Aids Incubation Times and Reporting Lags.  _Statistica Sinica_ , 1(1):19–32, 1991.
  * Keiding and Moeschberger (1992) Niels Keiding and Melvin Moeschberger.  Independent Delayed Entry.  In _Survival Analysis: State of the Art_ , pages 309–326. Springer Netherlands, 1992.  10.1007/978-94-015-7983-4_18.
  * Lawless (1994) J. F. Lawless.  Adjustments for Reporting Delays and the Prediction of Occurred but Not Reported Events.  _The Canadian Journal of Statistics / La Revue Canadienne de Statistique_ , 22(1):15, 1994.  ISSN 03195724.  10.2307/3315820.
  * Midthune et al. (2005) Douglas N. Midthune, Michael P. Fay, Limin X. Clegg, and Eric J. Feuer.  Modeling reporting delays and reporting corrections in cancer registry data.  _Journal of the American Statistical Association_ , 100(469):61–70, 2005.  ISSN 01621459.  10.1198/016214504000001899.
  * Morgan and Curran (1986) W. M. Morgan and J. W. Curran.  Acquired immunodeficiency syndrome: Current and future trends.  _Public Health Reports_ , 101(5):459–465, 1986\.  ISSN 00333549.
  * Noufaily et al. (2015) Angela Noufaily, Yonas Ghebremichael-Weldeselassie, Doyo Gragn Enki, Paul Garthwaite, Nick Andrews, André Charlett, and Paddy Farrington.  Modelling reporting delays for outbreak detection in infectious disease data.  _Journal of the Royal Statistical Society. Series A: Statistics in Society_ , 178(1):205–222, 2015.  ISSN 1467985X.  10.1111/rssa.12055.
  * Noufaily et al. (2016) Angela Noufaily, Paddy Farrington, Paul Garthwaite, Doyo Gragn Enki, Nick Andrews, and Andre Charlett.  Detection of Infectious Disease Outbreaks From Laboratory Data With Reporting Delays.  _Journal of the American Statistical Association_ , 111(514):488–499, 2016.  ISSN 1537274X.  10.1080/01621459.2015.1119047.
  * Rosenberg (1990) Philip S. Rosenberg.  A simple correction of AIDS surveillance data for reporting delays.  _Journal of Acquired Immune Deficiency Syndromes_ , 3(1):49–54, 1990.  ISSN 10779450.
  * Wang (1992) Mei-Cheng Wang.  The Analysis of Retrospectively Ascertained Data in the Presence of Reporting Delays.  _Journal of the American Statistical Association_ , 87(418):397, 1992.  ISSN 01621459.  10.2307/2290270.
  * Weinberger et al. (2020) Daniel M. Weinberger, Jenny Chen, Ted Cohen, Forrest W. Crawford, Farzad Mostashari, Don Olson, Virginia E. Pitzer, Nicholas G. Reich, Marcus Russi, Lone Simonsen, Anne Watkins, and Cecile Viboud.  Estimation of Excess Deaths Associated with the COVID-19 Pandemic in the United States, March to May 2020.  _JAMA Internal Medicine_ , 06520(May):E1–E9, 2020\.  ISSN 21686114.  10.1001/jamainternmed.2020.3391.  URL<EMAIL_ADDRESS>
  * White et al. (2009) Laura Forsberg White, Jacco Wallinga, Lyn Finelli, Carrie Reed, Steven Riley, Marc Lipsitch, and Marcello Pagano.  Estimation of the reproductive number and the serial interval in early phase of the 2009 influenza A/H1N1 pandemic in the USA.  _Influenza and other Respiratory Viruses_ , 3(6):267–276, nov 2009.  ISSN 17502640.  10.1111/j.1750-2659.2009.00106.x.
  * Zhao and Zhou (2010) Xiao Bing Zhao and Xian Zhou.  Applying copula models to individual claim loss reserving methods.  _Insurance: Mathematics and Economics_ , 46(2):290–299, 2010.  ISSN 01676687.  10.1016/j.insmatheco.2009.11.001.
  * Zhao et al. (2009) Xiao Bing Zhao, Xian Zhou, and Jing Long Wang.  Semiparametric model for prediction of individual claim loss reserving.  _Insurance: Mathematics and Economics_ , 45(1):1–8, 2009.  ISSN 01676687.  10.1016/j.insmatheco.2009.02.009.