# Bayesian analysis of diffusion-driven multi-type epidemic models with
application to COVID-19

Lampros Bouranis∗,1, Nikolaos Demiris1,
Konstantinos Kalogeropoulos2, Ioannis Ntzoufras1

###### Abstract

We consider a flexible Bayesian evidence synthesis approach to model the age-
specific transmission dynamics of COVID-19 based on daily mortality counts.
The temporal evolution of transmission rates in populations containing
multiple types of individual is reconstructed via an appropriate dimension-
reduction formulation driven by independent diffusion processes. A suitably
tailored compartmental model is used to learn the latent counts of infection,
accounting for fluctuations in transmission influenced by public health
interventions and changes in human behaviour. The model is fitted to freely
available COVID-19 data sources from the UK, Greece and Austria and validated
using a large-scale seroprevalence survey in England. In particular, we
demonstrate how model expansion can facilitate evidence reconciliation at a
latent level. The code implementing this work is made freely available via the
Bernadette R package.

Keywords: Bayesian evidence synthesis; Brownian motion; COVID-19; Hamiltonian
Monte Carlo; Population epidemic model; Time-varying parameters.

##  1 Introduction

The COVID-19 outbreak caused by the Severe acute respiratory syndrome
coronavirus 2 (SARS-CoV-2), has led to developments in Bayesian infectious
disease modeling, allowing modelers to assess the impact of mitigation
strategies on transmission and quantify the burden of the pandemic (Flaxman et
al.,, 2020; Monod et al.,, 2021). The highly transmissible nature of COVID-19
has proven to be a challenge for health systems worldwide, particularly in the
acute phase of the pandemic, pre-vaccination. A number of large-scale
seroprevalence studies aiming to estimate the actual number of infections have
found severe under-ascertainment (Ward et al.,, 2021), which has been varying
in time and across countries. The level of under-ascertainment depends on the
national testing and tracing policies, the testing capacities and the impact
of false positives under different regimes. Especially at the early stages of
the pandemic, those tested were typically more likely to have been
hospitalised or were at higher risk of infection or death, leading to only a
proportion of infections being detected and reported (Li et al.,, 2020).
Consequently, methods that rely on reported counts of infections are expected
to yield biases in the inferred rates of transmission.

This work is driven by the challenges of under-ascertainment of COVID-19
infections and the presence of heterogeneity in type, relevance, and
granularity of the data. Building upon the approach by Flaxman et al., (2020)
and its extension to multiple age groups by Monod et al., (2021), we propose a
Bayesian evidence synthesis approach for the analysis of COVID-19 data, with
the aim to: (i) infer the true number of infections using age-stratified daily
COVID-19 attributable mortality counts; (ii) learn the age-specific
transmission rates; (iii) reconstruct the epidemic drivers from publicly
available data sources.

Chatzilena et al., (2022) extend the work of Dureau et al., (2013) allowing
for indirectly observed infections. We extend this framework by developing an
age-stratified transmission model with distinct contact rates between age
groups accounting for the presence of social structures. In particular, we
target the transmission rate matrix process, the dimension of which increases
quadratically with the number of age groups, and offer a dimension-reduction
formulation projecting to latent diffusion processes. This formulation
possesses several desirable characteristics in capturing the multiple group
disease transmission mechanism: (i) it is based on a natural decomposition
into the underlying biological and social components; (ii) it allows for
further evidence synthesis utilising information from contact surveys; (iii)
it is driven by a potentially non-scalar diffusion process to adequately
capture the temporal evolution of the transmission rates, as well as extrinsic
environmental factors such as non-pharmaceutical interventions (NPIs) and
climatic changes (Cazelles et al.,, 2021; Chatzilena et al.,, 2022; Gosh et
al.,, 2023); (iv) it facilitates model determination at a latent level,
performed here by appropriate model expansion.

We develop a Bayesian hierarchical method that integrates multiple data
sources which are publicly available across countries: age-stratified daily
COVID-19 attributable mortality counts and daily laboratory-confirmed COVID-19
infection counts, contact surveys and the age distribution of the population;
the contact survey data are used to delineate potential identifiability issues
(Britton,, 1998) at the unobserved infection rate level when those rates are
decomposed in their social and biological component. The uncertainty in
contact structure is expressed via suitable prior distributions. Our model
allows for borrowing information across age groups and time. We propose to do
this by representing the latent effective contact rate as independent
diffusion processes thus facilitating model realism and explainability.

The outline of the paper is as follows. Section 2 provides an overview of the
multiple data sources and a presentation of the components of the Bayesian
hierarchical model. In Section 3 we present the results of the empirical
analysis for Greece and discuss the effect of model expansion on resolving
prior-data conflict at a latent level. The proposed model is validated using
estimates from a large-scale seroprevalence survey in England. Results for
Austria are presented in the Supplementary Material. We conclude the paper in
Section 4 with final remarks.

##  2 Materials and Methods

###  2.1 Data sources

Our model is driven by the characteristics of SARS-CoV-2 transmission and
applied to observations of the epidemic in Greece, Austria and England. For
each of these countries, we collected three types of data related to SARS-
CoV-2 spread: (i) the age distribution of reported deaths; (ii) the age
distribution of all reported infections; (iii) the age distribution of the
general population for each country, $\mathbb{N}$. The study period for
Austria ranges from 2020-08-31 to 2021-03-28 (30 weeks). The study period for
Greece ranges from 2020-09-01 to 2021-03-29 (30 weeks). For model validation,
we analyse the period 2020-03-02 to 2021-09-27 (30 weeks) in England.

We adopted the country-specific synthetic contact matrices of Prem et al.,
(2021), which have been constructed upon adjusted contact patterns, based on
national demographic and socioeconomic characteristics. In the absence of a
synthetic contact matrix for England we used the one created for the United
Kingdom, a reasonable assumption due to the population overlap. The elements
of these matrices represent the daily average numbers of contacts between age
groups.

The severity of the disease is incorporated in the model via the age-
stratified infection fatality rate (IFR), informed by the REal-time Assessment
of Community Transmission-2 (REACT-2) national seroprevalence study of over
100,000 adults in England (Ward et al.,, 2021). Detailed references for each
data source are presented in Section S1 of the Supplementary Material.

###  2.2 Modeling framework

####  2.2.1 Evidence synthesis

The aforementioned data streams and expert knowledge were integrated into a
coherent modeling framework via a Bayesian evidence synthesis approach (De
Angelis and Presanis,, 2018). Following the works of Flaxman et al., (2020),
Monod et al., (2021) and Chatzilena et al., (2022), the modeling process is
separated into a latent epidemic process and an observation process in an
effort to reduce sensitivity to observation noise and to allow for more
flexibility in modeling different forms of data. Figure 1 shows the functional
relationships between data sources, modeled outputs and parameters in this
study. The components of the Bayesian hierarchical model are laid out in the
remainder of Section 2.2.

####  2.2.2 Diffusion-driven multi-type transmission process

The transmission of COVID-19 is modeled through an age-stratified
deterministic Susceptible-Exposed-Infectious-Removed (SEIR) compartmental
model (Anderson and May,, 1992) with Erlang-distributed latent and infectious
periods. In particular, we introduce Erlang-distributed stage durations in the
Exposed and Infected compartments and relax the mathematically convenient but
unrealistic assumption of exponential stage durations (Anderson and Watson,,
1980; Lloyd,, 2001) by assuming that each of the Exposed and Infected
compartments are defined by two stages, with the same rate of loss of latency
($\tau$) and infectiousness ($\gamma$) in both stages. Removed individuals are
assumed to be immune to reinfection for at least the duration of the study
period.

The population is stratified into $\alpha\in\\{1,\ldots,A\\}$ age groups and
the total size of the age group is denoted by
$\mathbb{N}_{\alpha}=S_{t}^{\alpha}+E_{1,t}^{\alpha}+E_{2,t}^{\alpha}+I_{1,t}^{\alpha}+I_{2,t}^{\alpha}+R_{t}^{\alpha}$,
where $S_{t}^{\alpha}$ represents the number of susceptible,
$E_{t}^{\alpha}=\sum_{j=1}^{2}E_{j,t}^{\alpha}$ represents the number of
exposed but not yet infectious,
$I_{t}^{\alpha}=\sum_{j=1}^{2}I_{j,t}^{\alpha}$ is the number of infected and
$R_{t}^{\alpha}$ is the number of removed individuals at time $t$ at age group
$\alpha$. The number of individuals in each compartment is scaled by the total
population $\mathbb{N}=\sum_{\alpha=1}^{A}\mathbb{N}_{\alpha}$, so that the
sum of all compartments equals to one (Grinsztajn et al.,, 2021). The latent
epidemic process is expressed by the following non-linear system of ordinary
differential equations (ODEs) (Kermack and McKendrick,, 1927)

$\begin{cases}\frac{dS^{\alpha}_{t}}{dt}&=-\lambda_{\alpha}(t)S^{\alpha}_{t},\\\
\frac{dE^{\alpha}_{1,t}}{dt}&=\lambda_{\alpha}(t)S^{\alpha}_{t}-\tau
E^{\alpha}_{1,t},\\\
\frac{dE^{\alpha}_{2,t}}{dt}&=\tau\left(E^{\alpha}_{1,t}-E^{\alpha}_{2,t}\right),\\\
\frac{dI^{\alpha}_{1,t}}{dt}&=\tau E^{\alpha}_{2,t}-\gamma
I^{\alpha}_{1,t},\\\
\frac{dI^{\alpha}_{2,t}}{dt}&=\gamma\left(I^{\alpha}_{1,t}-I^{\alpha}_{2,t}\right),\\\
\frac{dR^{\alpha}_{t}}{dt}&=\gamma I^{\alpha}_{2,t}\end{cases}$ (1)

where the mean latent and infectious periods are $d_{E}=\frac{2}{\tau}$,
$d_{I}=\frac{2}{\gamma}$, respectively. The number of new infections in age
group $\alpha$ at day $t$ is

$\Delta^{\text{infec}}_{t,\alpha}=\int_{t-1}^{t}\tau E^{\alpha}_{2,s}ds.$ (2)

The time-dependent force of infection $\lambda_{\alpha}(t)$ for age group
$\alpha\in\\{1,\ldots,A\\}$ is expressed as

$\lambda_{\alpha}(t)=\sum_{\alpha^{\prime}=1}^{A}\left[m_{\alpha,\alpha^{\prime}}(t)\frac{\left(I^{\alpha^{\prime}}_{1,t}+I^{\alpha^{\prime}}_{2,t}\right)}{\mathbb{N}_{\alpha^{\prime}}}\right],$

which is a function of the proportion of infectious individuals in each age
group $\alpha^{\prime}\in\\{1,\ldots,A\\}$, via the compartments
$I^{\alpha^{\prime}}_{1,t}$, $I^{\alpha^{\prime}}_{2,t}$ divided by the total
size of the age group $\mathbb{N}_{\alpha^{\prime}}$, and the time-varying
person-to-person transmission rate from group $\alpha$ to group
$\alpha^{\prime}$, $m_{\alpha,\alpha^{\prime}}(t)$. We parameterize the
transmission rate between different age groups
$(\alpha,\alpha^{\prime})\in\\{1,\ldots,A\\}^{2}$ by

$m_{\alpha,\alpha^{\prime}}(t)=\beta^{\alpha\alpha^{\prime}}_{t}\cdot
C_{\alpha,\alpha^{\prime}},$ (3)

breaking down the transmission rate matrix into its biological and social
components (Baguelin et al.,, 2013; Knock et al.,, 2021): the social component
is represented by the average number of contacts between individuals of age
group $\alpha$ and age group $\alpha^{\prime}$ via the contact matrix element
$C_{\alpha,\alpha^{\prime}}$; $\beta^{\alpha\alpha^{\prime}}_{t}$ is the time-
varying transmissibility of the virus, the probability that a contact between
an infectious person in age group $\alpha$ and a susceptible person in age
group $\alpha^{\prime}$ leads to transmission at time $t$.

The formulation below may be viewed as a stochastic extension to the
deterministic multi-type SEIR model, using diffusion processes for the
coefficients $\beta^{\alpha\alpha^{\prime}}_{t}$ in (3), driven by independent
Brownian motions

$\begin{cases}\beta^{\alpha\alpha^{\prime}}_{t}&=\exp(x^{\alpha\alpha^{\prime}}_{t}),\\\
x^{\alpha\alpha^{\prime}}_{t}\mid
x^{\alpha\alpha^{\prime}}_{t-1},\sigma^{\alpha\alpha^{\prime}}_{x}&\sim\operatorname{N}(x^{\alpha\alpha^{\prime}}_{t-1},(\sigma^{\alpha\alpha^{\prime}})^{2}_{x}),\\\
dx^{\alpha\alpha^{\prime}}_{t}&=\sigma^{\alpha\alpha^{\prime}}_{x}dW^{\alpha\alpha^{\prime}}_{t},\\\
dW^{\alpha\alpha^{\prime}}_{t}&\sim\operatorname{N}(0,d_{t}),\end{cases}$ (4)

with volatilities $\sigma^{\alpha\alpha^{\prime}}_{x}$, corresponding to the
case of little information on the shape of
$\beta^{\alpha\alpha^{\prime}}_{t}$. The volatility
$\sigma^{\alpha\alpha^{\prime}}_{x}$ plays the role of the regularizing
factor: higher values of $\sigma^{\alpha\alpha^{\prime}}_{x}$ lead to greater
changes in $\beta^{\alpha\alpha^{\prime}}_{t}$. The exponential transformation
avoids negative values which have no biological meaning (Dureau et al.,, 2013;
Cazelles et al.,, 2021). A major advantage of considering a diffusion process
for modeling $\beta^{\alpha\alpha^{\prime}}_{t}$ is its ability to capture and
quantify the randomness of the underlying transmission dynamics, which is
particularly useful when the dynamics are not completely understood. The
diffusion process accounts for fluctuations in transmission that are
influenced by non-modeled phenomena, such as new variants, mask-wearing
propensity, etc. The diffusion process also allows for capturing the effect of
unknown extrinsic factors on the age-stratified force of infection, for
monitoring of the temporal evolution of the age-specific transmission rate
without the implicit inclusion of external variables and for tackling non-
stationarity in the data.

We propose a diffusion-driven multi-type latent transmission model which
assigns independent Brownian motions to
$\log(\beta^{11}_{t}),\log(\beta^{22}_{t}),\ldots,\log(\beta^{AA}_{t})$ with
respective age-stratified volatility parameters
$\sigma_{x,\alpha},\alpha\in\\{1,\ldots,A\\}$, for reasons of parsimony and
interpretability. The contact matrix is scaled by the age-stratified
transmissibility in order to obtain the transmission rate matrix process

$m^{\text{MBM}}_{\alpha,\alpha^{\prime}}(t)=\beta^{\alpha\alpha^{\prime}}_{t}\cdot
C_{\alpha,\alpha^{\prime}}\equiv\beta^{\alpha\alpha}_{t}\cdot
C_{\alpha,\alpha^{\prime}},$ (5)

under the assumption
$\beta^{\alpha\alpha^{\prime}}_{t}\equiv\beta^{\alpha\alpha}_{t},\alpha\neq\alpha^{\prime}$.
Henceforth, we term the model in (5) by MBM (standing for Multi Brownian
Motion). The MBM model can also be viewed as factor analysis, offering
dimension reduction from all the elements of the transmission rate matrix
process to the Brownian motions $\beta_{t}^{\alpha\alpha}$ for
$\alpha\in\\{1,\ldots,A\\}$. An appealing feature of this factor analysis
representation is that it separates the contact matrix which reflects social
components of the model and can be informed from survey data. We provide
details on the formulation of the contact matrix in the Supplementary
Material.

Further reduction in the dimension of the transmission rate matrix process can
be achieved by the formulation

$m^{\text{SBM}}_{\alpha,\alpha^{\prime}}(t)=\beta^{\alpha\alpha^{\prime}}_{t}\cdot
C_{\alpha,\alpha^{\prime}}\equiv\beta_{t}\cdot C_{\alpha,\alpha^{\prime}},$
(6)

under the assumption $\beta^{\alpha\alpha^{\prime}}_{t}\equiv\beta_{t}$, where
$\beta_{t}$ is the time-varying transmissibility of the virus, the probability
that a contact between an infectious person and a susceptible person leads to
transmission at time $t$. Henceforth, we term the model in (6) by SBM
(standing for Single Brownian Motion), which represents a nested model to MBM.

In related studies, Birrell et al., (2021) estimated the transmission rate
between different age groups considering time-dependent (weekly-varying)
contact matrices, not allowing the age-stratified transmissibility to vary
unconstrained in time. Knock et al., (2021) considered a time-independent
contact matrix; the transmission rate between different age groups was
parameterised as in (6), however the virus transmissibility was assumed to be
piecewise linear with multiple change points corresponding to major
announcements and changes in COVID-19 related policy.

####  2.2.3 Observation process

Denote the number of observed deaths on day $t=1,\ldots,T$ in age group
$\alpha\in\\{1,\ldots,A\\}$ by $y_{t,\alpha}$. A given infection may lead to
observation events (i.e deaths) in the future. A link between $y_{t,\alpha}$
and the expected number of new age-stratified infections is established via
the function

$d_{t,\alpha}=\mathbb{E}[y_{t,\alpha}]=\widehat{\text{IFR}}_{\alpha}\times\sum_{s=1}^{t-1}h_{t-s}\Delta^{\text{infec}}_{s,\alpha}$

on the new expected age-stratified mortality counts, $d_{t,\alpha}$, the
estimated age-stratified infection fatality rate,
$\widehat{\text{IFR}}_{\alpha}$, and the infection-to-death distribution $h$,
where $h_{s}$ gives the probability the death occurs on the $s^{th}$ day after
infection, is assumed to be gamma distributed with mean 24.2 days and
coefficient of variation 0.39 (Flaxman et al.,, 2020; Monod et al.,, 2021),
that is

$h\sim\operatorname{Gamma}(6.29,0.26).$ (7)

We allow for over-dispersion in the observation processes to account for noise
in the underlying data streams, for example due to day-of-week effects on data
collection (Stoner et al.,, 2020; Knock et al.,, 2021), and link
$d_{t,\alpha}$ to $y_{t,\alpha}$ through an over-dispersed count model (Hauser
et al.,, 2020; Birrell et al.,, 2021; Seaman et al.,, 2022)

$y_{t,\alpha}\mid\theta\sim\operatorname{NegBin}\left(d_{t,\alpha},\xi_{t,\alpha}\right),$
(8)

where $\xi_{t,\alpha}=\frac{d_{t,\alpha}}{\phi}$, such that
$\mathbb{V}[y_{t,\alpha}]=d_{t,\alpha}(1+\phi)$. The log-likelihood of the
observed deaths is given by

$\ell^{Deaths}(y\mid\theta)=\sum_{t=1}^{T}\sum_{\alpha=1}^{A}\text{logNegBin}\left(y_{t,\alpha}\mid
d_{t,\alpha},\xi_{t,\alpha}\right),$

where $y\in\mathbb{R}^{T\times A}_{0,+}$ are the surveillance data on deaths
for all time-points and age groups and the parameter vector corresponds to
either $\theta^{SBM}$ or $\theta^{MBM}$, defined in Section 2.3. The
observation process in this work does not additionally model the age-
stratified infection counts, for reasons discussed in Section 1.

###  2.3 Parameter estimation

Bayesian estimation of the parameters of the SBM and MBM models was performed
using Markov chain Monte Carlo (MCMC). Computational details are provided in
the Supplementary Material. In the following we discuss the challenging
aspects of parameter estimation and present inference approximations used in
the implementation in order to improve mixing and computational efficiency.

Estimation of the new daily infection counts in (2) requires the solution of
the non-linear system of ODEs in (1) coupled with the stochastic differential
equation in (4), which together are cast as a hypo-elliptic diffusion which is
intractable. We address the intractability of the hypo-elliptic diffusion by
adopting the data augmentation framework of Dureau et al., (2013) as a means
to infer the latent sample path $x$ of the diffusion which is an infinite-
dimensional object (Kalogeropoulos et al.,, 2010) and is indirectly observed
through the time evolution of the disease states. The task of inference for
the transmission rate(s) via diffusions is particularly challenging; modeling
its time evolution at a more granular time-scale (i.e. at the time-points of
the observations) would naturally increase the cost of the computational
effort; to reduce the dimensionality of the stochastic process we split the
study period into $k=1,\ldots,K$ weeks and denote by $k_{t}$ the week that
time point $t$ falls into. We assumed that the transmissibility of the virus,
$\beta_{k_{t}}$, remains constant between subsequent weeks $[k_{t},k_{t}+1)$
and employed time-discretization via an Euler-Maruyama approximation of the
latent sample path $x$. Switching the notation from $t$ to $k_{t}$, this
implies $x_{k_{t}}\mid
x_{k_{t}-1},\sigma_{x}\sim\operatorname{N}(x_{k_{t}-1},\sigma^{2}_{x})$ for
the SBM model, which may be viewed as a prior distribution on $x_{k_{t}}$ and
decreases the dimensionality of the stochastic process to $K+1<T+1$.

Another challenging aspect of parameter estimation is the estimation of the
volatility $\sigma_{x}$ which is a top-level parameter in the Bayesian
hierarchical model; due to the multiple levels of hierarchy in the evidence
synthesis model (Figure 1) and the lack of information about it, this
difficulty is reflected in the lower effective sample sizes in all of our
analyses. Subsequently, this increases the difficulty to break down the
transmission rate in (3) into its biological and social components. Finally,
for this multivariate ODE time series model, the solution to (1) is
approximated numerically via the Trapezoidal rule. Higher-order numerical ODE
solvers can be considered, at the cost of increased computational effort.

The Bayesian hierarchical model is implemented via a dynamic Hamiltonian Monte
Carlo algorithm (Betancourt,, 2018) which utilizes gradient information to
efficiently explore high-dimensional parameter spaces or complex posterior
geometries and to obtain a sample from the posterior distribution of the model
parameters given the observed data. The size of the parameter space increases
with longer time horizons due to the involvement of the diffusion process in
(4) at weekly level, making such a gradient-based MCMC method particularly
suitable for the implementation of the Bayesian hierarchical model proposed in
the work; simultaneously, this poses a challenge for ODE models, for which
computing the gradient is computationally significantly more demanding and
subtle than the plain likelihood evaluation (Timonen et al.,, 2022). The
increasing computational cost as a consequence of increasing volumes of data
for ODE models is a point raised by Birrell et al., (2021) and Gosh et al.,
(2023) and highlights the need for further developments in the area.

##  3 Data Analysis

###  3.1 The COVID-19 pandemic data

The proposed methodology is illustrated on freely available data from the
COVID-19 pandemic in Greece between September 2020 and March 2021, where
$K=30$ weeks. The population was divided in three age groups,
$\\{0-39,40-64,65+\\}$. During the study period, a national lockdown was
implemented in November 2020 until the end of January, 2021. Some of the NPIs
were relaxed at New Year’s Day (Table S3 of the Supplementary Material). In
this section we focus on the SBM model which accounts for the inherent
variability of contacts in the population as a result of the response of
individuals to NPIs which were taken to reduce transmission, expressing the
uncertainty in contact structure via prior distributions on the elements of
the random contact matrix $C$, which is constrained to not vary in time.

However, the assumption made in (6) restricts the age-stratified effective
contact rates to be expressed by an overall effective contact rate for the
population, which is dominated by the effective contact rate of the age group
that drives the transmission of the virus, at the time-points of the
observations (Figure 2, panel A). This non-realistic assumption does not
properly account for the change in the behaviour of individuals at each age
group over time as a result of the implementation of NPIs during the study
period. As previously discussed in Section 2.3, the lack of sufficient
information about the volatility $\sigma_{x}$ increases the difficulty to
break down the transmission rate in (6) into its biological and social
components, creating an identifiability issue. Therefore, any differences in
the virus transmission procedure across strata are attributed to the random
contact matrix elements $C_{\alpha,\alpha^{\prime}}$, justifying the notable
differences between the prior and the posterior distribution of the contact
matrix (Figure 2, panel B). Consequently, the age-stratified transmission rate
trajectories (Figure 2, panel C) only differ in terms of magnitude at the
time-points of the observations.

The analysis of Austrian mortality counts revealed the same view regarding the
ability of the SBM model to capture the age-stratified trends in SARS-CoV-2
transmission (Figure S5 of the Supplementary Material). These findings suggest
that the SBM model is not flexible enough to accommodate for age-stratified
trends in SARS-CoV-2 transmission, motivating the need for a more flexible
model which can better estimate the volatility of the diffusion process.

###  3.2 Model expansion

The assessment of the fidelity of the SBM model to the data at a latent level
revealed an issue which may be seen as prior-data conflict at a latent level.
We resolve this issue by expanding the SBM model to the MBM model in the
spirit of Gustafson, (2005) and inspecting the effect on the contact survey
data. This type of evidence reconciliation is not trivial to assess and it
appears that inspecting the prior to posterior updates for the random but
time-independent contact matrix yields a reasonable approach to surmount the
conflict.

The mortality counts that are available from the national surveillance systems
consistently indicate an increasing number of COVID-19 attributable deaths
with age; for the case of Greece the respective time series for the younger
age group, i.e. $\\{0-39\\}$, is sparsely-informed, while the elder age groups
$\\{40-64,65+\\}$ are well-informed; a similar observation can be made for the
Austrian dataset (Figures S1-S4 of the Supplementary Material). In the context
of diffusion-driven transmission models it would be more sensible to group
younger ages together so as to avoid sparsely-informed time series which are
more difficult to fit and instead allow for more groupings in elder ages, i.e.
by 5- or 10-year intervals, depending on the information that becomes
available from the national surveillance systems. In such a case, where the
number of age groups $A$ is small, the mortality data do not provide much
information on the top-level volatilities
$\sigma_{x,\alpha},\alpha\in\\{1,\ldots,A\\}$ of the MBM Bayesian hierarchical
model, reflecting the difficulty in the estimation of the volatility
$\sigma_{x}$ involved in the SBM model. Similarly to the SBM model, posterior
inference of the volatilities $\sigma_{x,\alpha},\alpha\in\\{1,\ldots,A\\}$
has resulted in lower effective sample sizes compared to the rest of the MBM
model parameters.

The assumption of independent diffusions for the biological components of the
transmission rate in (5) allows for better reconstructing of the age-
stratified drivers of transmission (Figure 3, panels A and C). The order of
magnitude of the transmission rates changes between the SBM and MBM models
because of the change in the order of magnitude in the posterior number of
contacts (Figures 2 and 3, panel B). Overall, the expanded transmission model
offers better interpetability and more flexibility than the SBM model in
inferring dynamic transmission rates and enables the reconstruction of the
age-stratified drivers of transmission. We discuss a quantitative comparison
of the two models using information criteria in Section S3 of the
Supplementary Material.

The MBM model is validated in the period 2020-03-02 to 2021-09-27 ($K=30$
weeks) using the estimated age-stratified numbers of cumulative infections in
England from the REACT-2 seroprevalence survey (Ward et al.,, 2021). In the
analysis the population was divided in three age groups,
$\\{0-39,40-59,60+\\}$. REACT-2 quantified the total number of infected
individuals in England and the level of under-ascertainment of infections in
the population up to mid-July 2020; it was selected for model validation based
on its careful design, large representative sample size and timing. The
estimated counts of cumulative infections from REACT-2 were adjusted by the
age distribution of the population based on the three age groups
$\\{0-39,40-59,60+\\}$. Figure 4 shows the agreement of the model estimates
with the estimates from REACT-2 and demonstrates the ability of the MBM model
to capture the level of under-ascertainment of infections over time and by age
group.

###  3.3 Effectiveness of interventions

We assess the interventions in transmission (Table S4 of the Supplementary
Material) for Greece based on the findings of the MBM transmission model. The
findings for Austria are discussed in the Supplementary Material.

The trajectory of the posterior transmission rate for individuals aged $60+$
(Figure 3, panel C) resembles the respective trajectory shown in Figure 2
Panel C, thus stressing the role of this group to the evolution of the
epidemic. The transmission rate of the virus reached its highest levels in the
third week of October 2020; the largest magnitude is demonstrated for the age
group $\\{60+\\}$ (Figure 3, panel C), providing evidence to support public
health authorities in taking appropriate measures to decrease transmission
within and between-groups. The consecutive decline was steeper for individuals
aged $40+$ in the following two weeks. Interventions to limit all
indoor/outdoor mass/public gatherings until the end of January 2021 were
particularly effective in individuals aged $40+$, keeping virus transmission
to low levels, at about 10%. The interventions implemented at the end of
January 2021 to limit outdoor mass/public gatherings of over 100 participants
were successful in controlling the gradual increase in transmission in
individuals aged $40+$. Partial relaxation of the measure regarding the
closure of day cares, primary and secondary schools in early February 2021 for
a month appears to have contributed in controlling transmission.

##  4 Discussion

In this paper we have presented an epidemic model driven by the distinct
features of COVID-19 data. We proposed a flexible Bayesian evidence synthesis
framework that enables a data-driven approach to inferring the mechanism
governing COVID-19 transmission, based on diffusion processes that are
a-priori independent. The modularity of the Bayesian approach allows for
assessing fidelity to the data at a latent level and resolving the
corresponding prior-data conflict via expanding the model to incorporate
distinct diffusions for each age group.

There are several advantages associated with our approach. Our model is
primarily driven by the reported daily COVID-19 attributable mortality counts;
a key strength of our Bayesian evidence synthesis framework is that multiple
data sources which are publicly available across countries are integrated to
provide a robust overall picture of the epidemic nationwide: daily COVID-19
attributable mortality counts, daily laboratory-confirmed COVID-19 infection
counts and the age distribution of the population. This increases the
applicability of the model to multiple countries, since these data sources are
made available by the national surveillance systems and we do not rely on
bespoke surveys. We develop data-driven estimation of the transmissibility of
the virus via mutually independent diffusion processes thus removing the need
for additional information on the timing of specific interventions and
hypotheses about their impact on transmission. Instead, more general
assumptions need to be made, such as on the smoothness of the latent sample
path of the age-stratified diffusion that we wish to infer.

Our approach allows policymakers to assess the effect of NPIs on each age
group. The model simultaneously estimates multiple metrics of outbreak
progression: by modeling mortality counts instead of inferring only from
confirmed infections, the Bayesian hierarchical model is more robust to
changes in COVID-19 testing policies and the ability of national healthcare
systems to detect COVID-19 infections, i.e. through contact tracing, and
provides valuable information about the age-stratified transmission rates,
latent number of infections, reporting ratio and the effective reproduction
number (Sections S4-S7 of the Supplementary Material), where uncertainty is
accounted for naturally via MCMC.

Our modeling approach is subject to limitations. Age-stratified counts of
vaccinations are not accounted for in the latent transmission epidemic model
in the form of covariates, but this does not pose an important limitation for
the considered time period, with low counts of administered vaccinations until
the end of March 2021. In its current form the proposed hierarchical model is
better suited to the pre-vaccine era. In addition, the model is limited by
variations in the reporting procedure of deaths and mortality definitions
across time and countries (Section S1 of the Supplementary Material); the
dates of daily reported deaths may deviate from the actual dates of deaths.
The use of diffusion processes can be viewed as an attempt to model the noise
within the observational model and mitigate such issues. Similar to the
Bayesian models of Monod et al., (2021) and Wistuba et al., (2022), we make
parametric assumptions regarding the distribution of the time between
infection and reported death in (7) (i.e. it is assumed to be the same across
age groups and constant during the study period). While factors like hospital
capacity and vaccine efficacy are expected to affect the infection-to-death
distribution, such parametric assumptions might not be generally transferable
to countries with different reporting characteristics and healthcare systems.

Some further limitations are commonplace in the literature. We model a closed
environment with no infections imported from outside of the population.
Additionally, our model does not explicitly account for hospital-acquired
infections which may have contributed to overall transmission and to
persistence of infection in periods of high transmission. Age-specific counts
of hospitalizations and intensive care unit cases (where publicly available)
could help relax the parametric assumptions regarding the infection-to-death
distribution and assist in the estimation of the course of the COVID-19
pandemic. The age-stratified IFR is a crucial component of the observation
model as it enables estimation of the expected number of new infections from
the observed mortality counts. An assumption that has been made is that the
age-stratified IFRs remain constant during the study period. However, these
are expected to vary in time and across populations due to a number of
factors, such as SARS-CoV-2 variants, vaccine efficacy, the age distribution
of the population, the age distribution of infections, underlying health
conditions, access to medication, the overall burden to the healthcare system
and other factors (Brazeau et al.,, 2022). Detailed such data could be
incorporated leading to a dynamic age-stratified IFR.

In future work we shall work on two directions to expand our model. The first
will aim to incorporate dynamic age-stratified IFR estimates and integrate
further healthcare surveillance data like age-stratified counts of
vaccinations, hospitalizations and intensive care unit cases where these data
are available. The second extension relates to multi-task Gaussian Processes
which have been recently implemented for a different class of Bayesian
infectious disease models (Seymour et al.,, 2022) could amend our model via
exchangeable diffusion processes, reflecting a-priori beliefs that there is a
shared structure between the dynamic transmission rates for individuals of
different age groups.

## Acknowledgments

The authors would like to thank Anastasia Chatzilena, Theodore Kypraios and
David Rossell for their helpful and constructive comments. Lampros Bouranis is
supported by the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie grant agreement No 101027218.
Conflict of Interest: None declared.

## Supplementary Material

The reader is referred to the online Supplementary Material for additional
information regarding this work. The Bernadette package (Bouranis,, 2023) for
R (R Core Team,, 2023) implements the methodology in this paper. The R code
and documentation to reproduce the analysis are available at
https://github.com/bernadette-eu/indepgbm.

## References

  * Anderson and Watson, (1980) Anderson, D. and Watson, R. (1980).  On the spread of a disease with gamma distributed latent and infectious periods.  Biometrika, 67:191–198.
  * Anderson and May, (1992) Anderson, R. and May, R. (1992).  Infectious diseases of humans: dynamics and control.  Oxford university press, Oxford, United Kingdom.
  * Baguelin et al., (2013) Baguelin, M., Flasche, S., Camacho, A., Demiris, N., Miller, E., and Edmunds, W. (2013).  Assessing optimal target populations for influenza vaccination programmes: An evidence synthesis and modelling study.  PLoS Med., 10(10):e1001527.
  * Betancourt, (2018) Betancourt, M. (2018).  A conceptual introduction to Hamiltonian Monte Carlo.  https://arxiv.org/abs/1701.02434.
  * Birrell et al., (2021) Birrell, P., Blake, J., van Leeuwen, E., Gent, N., and De Angelis, D. (2021).  Real-time nowcasting and forecasting of COVID-19 dynamics in England: the first wave.  Phil Trans R Soc B, 376:20200279.
  * Bouranis, (2023) Bouranis, L. (2023).  Bernadette: Bayesian Inference and Model Selection for Stochastic Epidemics.  R package version 1.1.4.
  * Brazeau et al., (2022) Brazeau, N., Verity, R., Jenks, S., Fu, H., Whittaker, C. andWinskill, P., Dorigatti, I., Walker, P., Riley, S., Schnekenberg, R., Hoeltgebaum, H., Mellan, T., Mishra, S., Unwin, H., Watson, O., Cucunuba, Z., Baguelin, M., Whittles, L., Bhatt, S., Ghani, A., Ferguson, N., and Okell, L. (2022).  Estimating the COVID-19 infection fatality ratio accounting for seroreversion using statistical modelling.  Commun. Med., 2(54).
  * Britton, (1998) Britton, T. (1998).  Estimation in multitype epidemics.  J. R. Statist. Soc. B, 60(4):663–679.
  * Cazelles et al., (2021) Cazelles, B., Champagne, C., Nguyen-Van-Yen, B., Comiskey, C., Vergu, E., and Roche, B. (2021).  A mechanistic and data-driven reconstruction of the time-varying reproduction number: Application to the COVID-19 epidemic.  PLos Comput Biol, 17(7):1–20.
  * Chatzilena et al., (2022) Chatzilena, A., Demiris, N., and Kalogeropoulos, K. (2022).  A modelling framework for the analysis of the transmission of SARS-CoV2.  https://arxiv.org/abs/2203.03773.
  * De Angelis and Presanis, (2018) De Angelis, D. and Presanis, A. (2018).  Analysing multiple epidemic data sources.  In Held, L., Hens, N., O’Neill, P., and Wallinga, J., editors, Handbook of Infectious Disease Data Analysis. Chapman & Hall/CRC, Florida, USA.
  * Dureau et al., (2013) Dureau, J., Kalogeropoulos, K., and Baguelin, M. (2013).  Capturing the time-varying drivers of an epidemic using stochastic dynamical systems.  Biostatistics, 14(3):541–555.
  * Flaxman et al., (2020) Flaxman, S., Mishra, S., Gandy, A., Unwin, H., Mellan, T., Coupland, H., Whittaker, C., Zhu, H., Berah, T., and Eaton, J. (2020).  Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe.  Nature, 584(7820):257–261.
  * Gosh et al., (2023) Gosh, S., Birrell, P., and De Angelis, D. (2023).  An approximate diffusion process for environmental stochasticity in infectious disease transmission modelling.  PLoS Comput Biol, 19(5):e1011088.
  * Grinsztajn et al., (2021) Grinsztajn, L., Semenova, E., Margossian, C., and Riou, J. (2021).  Bayesian workflow for disease transmission modeling in stan.  Statist Med, 40(27):6209–6234.
  * Gustafson, (2005) Gustafson, P. (2005).  On Model Expansion, Model Contraction, Identifiability and Prior Information: Two Illustrative Scenarios Involving Mismeasured Variables.  Stat Sci, 20(2):111–140.
  * Hauser et al., (2020) Hauser, A., Counotte, M., Margossian, C., Konstantinoudis, G., Low, N., Althaus, C., and Riou, J. (2020).  Estimation of SARS-CoV-2 mortality during the early stages of an epidemic: A modeling study in Hubei, China, and six regions in Europe.  PLoS Med, 17(7):1–17.
  * Kalogeropoulos et al., (2010) Kalogeropoulos, K., Roberts, G., and Dellaportas, P. (2010).  Inference for stochastic volatility models using time change transformations.  Ann. Statist., 38(2):784–807.
  * Kermack and McKendrick, (1927) Kermack, W. and McKendrick, A. (1927).  A contribution to the mathematical theory of epidemics.  Proc. R. Soc. Lond. A, 115:700–721.
  * Knock et al., (2021) Knock, E., Whittles, L., Lees, J., Perez-Guzman, P., Verity, R., FitzJohn, R., Gaythorpe, K., Imai, N., Hinsley, W., Okell, L., Rosello, A., Kantas, N., Walters, C., Bhatia, S., Watson, O., Whittaker, C., Cattarino, L., Boonyasiri, A., Djaafara, B., Fraser, K., Fu, H., Wang, H., Xi, X., Donnelly, C., Jauneikaite, E., Laydon, D., White, P., Ghani, A., Ferguson, N., Cori, A., and Baguelin, M. (2021).  Key epidemiological drivers and impact of interventions in the 2020 SARS-CoV-2 epidemic in England.  Sci. Transl. Med., 13(602):eabg4262.
  * Li et al., (2020) Li, R., Pei, S., Chen, B., Song, Y., Zhang, T., Yang, W., and Shaman., J. (2020).  Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (sars-cov-2).  Science, 368(6490):89–493.
  * Lloyd, (2001) Lloyd, A. (2001).  Realistic distributions of infectious periods in epidemic models: Changing patterns of persistence and dynamics.  Theor. Popul. Biol., 60:59–71.
  * Monod et al., (2021) Monod, M., Blenkinsop, A., Xi, X., Hebert, D., Bershan, S., Tietze, S., Bradley, V., Chen, Y., Coupland, H., Filippi, S., Ish-Horowicz, J., McManus, M., Mellan, T., Gandy, A., Hutchinson, M., Unwin, J., Vollmer, C., Michaela, A., Weber, S., Zhu, H., Bezancon, A., Ferguson, N., Mishra, S., Flaxman, S., Bhatt, S., and Ratmann, O. (2021).  Age groups that sustain resurging COVID-19 epidemics in the United States.  Science, 371(6536):eabe8372.
  * Prem et al., (2021) Prem, K., Zandvoort, K., Klepac, P., Eggo, R., Davies, N., Cook, A., and Jit, M. (2021).  Projecting contact matrices in 177 geographical regions: An update and comparison with empirical data for the COVID-19 era.  PLoS Comput Biol, 17(7):1–19.
  * R Core Team, (2023) R Core Team (2023).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Seaman et al., (2022) Seaman, S., Samartsidis, P., Kall, M., and De Angelis, D. (2022).  Nowcasting covid-19 deaths in england by age and region.  J. R. Statist. Soc. C, 71(5):1266–1281.
  * Seymour et al., (2022) Seymour, R., Kypraios, T., and O’Neill, P. (2022).  Bayesian nonparametric inference for heterogeneously mixing infectious disease models.  Proc. Natl. Acad. Sci., 119(10):e2118425119.
  * Stoner et al., (2020) Stoner, O., Economou, T., and Halliday, A. (2020).  A Powerful Modelling Framework for Nowcasting and Forecasting COVID-19 and Other Diseases.  https://arxiv.org/abs/1912.05965.
  * Timonen et al., (2022) Timonen, J., Siccha, N., Bales, B., Lähdesmäki, H., and Aki Vehtari, A. (2022).  An importance sampling approach for reliable and efficient inference in Bayesian ordinary differential equation models.  https://arxiv.org/abs/2205.09059.
  * Ward et al., (2021) Ward, H., Atchison, C., Whitaker, M., Ainslie, K., Elliott, J., Okell, L., Redd, R., Ashby, D., Donnelly, C., Barclay, W., Darzi, A., Cooke, G., Riley, S., and Elliott, P. (2021).  Sars-cov-2 antibody prevalence in england following the first peak of the pandemic.  Nat. Commun., 12:905.
  * Wistuba et al., (2022) Wistuba, T., Mayr, A., and Staerk, C. (2022).  Estimating the course of the COVID-19 pandemic in Germany via spline-based hierarchical modelling of death counts.  Sci. Rep., 12(9784).

Figure 1: Functional relationships between data sources (rectangles), modeled
outputs (ovals) and parameters (hexagons). Figure 2: Greece - Components of
the age-stratified transmission dynamics under the SBM transmission model.
Estimated median posterior trajectory (50% credible interval, CrI) trajectory
of the transmissibility of SARS-CoV-2 (panel A); prior and posterior
distributions of each element of the contact matrix (panel B); estimated
median posterior trajectory (50% credible interval, CrI) of the age-stratified
transmission rate (panel C).  Figure 3: Greece - Components of the age-
stratified transmission dynamics under the MBM transmission model. Estimated
median posterior trajectory (50% credible interval, CrI) of the
transmissibility of SARS-CoV-2 (panel A); prior and posterior distributions of
each element of the contact matrix (panel B); estimated median posterior
trajectory (50% credible interval, CrI) of the age-stratified transmission
rate (panel C).  Figure 4: External validation - England. Posterior median
age-stratified cumulative infections (50% & 95% credible intervals, CrI) under
the MBM model. The age-adjusted estimated counts and the respective 95%
confidence intervals of cumulative infections based on REACT-2 (Ward et al.,,
2021) are represented by horizontal lines. The vertical line corresponds to
middle of July 2020.


# Properties of $k$-Descending Trees

Agniv Sarkar and Eric Severson<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>973 Mission St, San Francisco, CA

(Date: 2021-2022)

###### Abstract.

This was research presented at the Worldwide Federation of National Math
Competitions in Bulgaria in 2022.

For any real-valued $k>1$, we consider the tree rooted at 0, where each
positive integer $n$ has parent $\lfloor\frac{n}{k}\rfloor$. The average
number of children per node is $k$, thus this definition gives a natural way
to extend $k$-ary trees to irrational $k$. We focus on the sequence $r_{d}$:
the count of nodes at depth $d$.

We first prove there exists some constant $\rho(k)$ such that
$r_{d}\sim\rho(k)\cdot k^{d}$. We then study a family of values
$k=\frac{a+\sqrt{a+4b}}{2}$, where we prove the sequence satisfies the exact
recurrence $r_{d}=a\cdot r_{d-1}+b\cdot r_{d-2}$. This generalizes a special
case when $k$ is the golden ratio and $r_{d}$ is the Fibonacci sequence.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
    1. 2.1 Notation
    2. 2.2 Definitions
  3. 3 Asymptotics of $r_{d}$
  4. 4 $k$-Trees with Closed Form $\rho(k)$
  5. 5 Conclusion
  6. A Tree Generation Code
  7. B $\rho$ Code
  8. C $\rho$ Raw Data
  9. D Grandparent Indicator Figures

##  1\. Introduction

This contains a piece of original research centered within graph and number
theory with a surprising connection to the Josephus Problem. The research was
carried out by myself (Agniv Sarkar) and my mentor (Eric Severson) throughout
the high school year of 2021-2022.

This research was done to observe patterns seen in the $\phi$-tree generated
with Definition 4. This became a very nice number theoretic problem, and when
we began to look at the asymptotics of these trees, we found that there was a
connection to the Josephus problem and calculating the solution to the
problem.

The trees themselves are most similar to a $k$-ary tree Definition 1.

######  Definition 1.

A $k$-ary tree is a rooted tree such that each node has no more than $k$
children. A complete $k$-ary tree is a rooted tree such that each node has
exactly $k$ children.

This tree is commonly used as a data structure in computer science, such as
through a Binary Search Tree, or a $2$-ary tree. However, the $k$ in the
definition does not generalize nicely to non integer $k$, and that is where
[2] defines a “rhythmic tree,” which is Definition 3.

######  Definition 2.

Let $p,q\in\mathbb{Z}$ such that $p>q\geq 1.$ Then,

  * •

Rhythm of directing parameter $(q,p)$ is a $q$-tuple of $r$ non-negative
integers whose sum is $p$.

$r=(r_{0},r_{1},\ldots,r_{q-1})\text{, and, }\sum_{i=0}^{q-1}r_{i}=p.$

  * •

A rhythm $r$ is valid if it also satisfies

$\forall k\in\\{0,1,\ldots,q-1\\},\sum_{i=0}^{j}r_{i}>j+1$

  * •

The growth rate of $r$ is the rational number $\frac{p}{q}$ or
$\frac{p^{\prime}}{q^{\prime}}$ where $p^{\prime}$ and $q^{\prime}$ are
quotients of $p$ and $q$ by their greatest common divisor, such that they are
coprime.

######  Definition 3.

Let $r=(r_{0},\ldots,r_{q-1})$ be a valid rhythm as given by Definition 2.
Then the rhythmic tree $\mathcal{I}_{r}$ generated by $r$ is defined by:

  * •

the root $0$ of $\mathcal{I}_{r}$ has $(r_{0}-1)$ children, which are the
notes $1,2,\ldots,$ and $(r_{0}-1)$.

  * •

for $n>0$, the node $n$ has $r_{n\mod q}$ children, which are the nodes
$(m+1),(m+2),\ldots,$ and $(m+r_{n\mod q})$ where $m$ is the largest child of
$(n-1)$.

In [2], they prove that if $r$ has rational but not integer growth rate, then
the paths to all vertices in the rhythmic tree $\mathcal{I}_{r}$ cannot be
verified with a finite automaton. If $r$ is a rhythm with integer growth rate,
then paths in the tree $\mathcal{I}_{r}$ can be described with a finite
automaton.

The reason that rhythmic trees are relevant and the previous statement about
their structure is that our new definition is a subset of all rhythmic trees
when $k$ is a rational number.

Also, [3] describes the Josephus problem, where we are given two numbers,
$n,q$. There are $n$ places arranged in a circle, and each $q$th person is
excused from the circle until only $1$ is left. The problem is to find out
what index the person who is left remains, which is simple for $q=2$, harder
for $q=3$, and there exists a recursive solution for a general $n,q$. This
paper defines the function $c(k).$

For some $k>1,$ fix the recursive sequence $f_{n}$ to be
$f_{0}=1,f_{n+1}=\lceil kf_{n}\rceil.$ Then, the function $c(k)$ is defined to
be the constant $c$ such that $f_{n}\sim c(k)k^{n}$ as $n$ trends to infinity.
We expand and prove more results about this function in Section 3.

##  2\. Preliminaries

###  2.1. Notation

$\mathbb{N}$ denotes the set of non-negative integers, $\mathbb{Q}$ the set of
rational numbers, and $\mathbb{R}$ the set of real numbers.
$\lfloor\cdot\rfloor$ and $\lceil\cdot\rceil$ denote the floor and ceiling
functions. For $x\in\mathbb{R}$, $\\{x\\}:=x-\lfloor x\rfloor$ denotes the
fractional part of $x$. Note that this means that
$\\{x\\}:\mathbb{R}\longrightarrow[0,1)$. $\sim$ denotes asymptotic
equivalence, where $f(n)\sim g(n)$ means
$\lim_{n\to\infty}\frac{f(n)}{g(n)}=1$.

###  2.2. Definitions

We first formally define a $k$-descending tree:

######  Definition 4.

For any $k\in(1,\infty)\subset\mathbb{R}$, the _$k$ -descending tree_, or
simply _$k$ -tree_ is the rooted tree with nodes in $\mathbb{N}$, where every
$n\in\mathbb{N}$ has the parent $\lfloor\frac{n}{k}\rfloor$, and $0$ is the
root node.

Fig. 1, Fig. 2, Fig. 3 were generated with Appendix A.

For the integer $k$ case, it generates a complete $k$-ary tree as shown in
Fig. 1.

Figure 1. $3$-descending tree — Used code to generate this tree.

Note that in Fig. 1 that it is the same as a rhythmic tree with rhythm $(3).$
This is true for all integers.

Then, for rational $k$, it generates another rhythmic tree as described in
[2]. The childcounts in Fig. 2 are periodic with pattern $2,1,\ldots$, meaning
that it is a rhythmic tree with rhythm $(2,1).$

Figure 2. $\dfrac{3}{2}$-descending tree — Coded it to generate this tree.

Then, when $k$ is irrational, the behavior becomes slightly more chaotic. This
specific example in Fig. 3 shows some particularly nice behavior due to our
choice of $k=\phi$.

Figure 3. $\frac{1+\sqrt{5}}{2}\approx 1.618$-descending tree — Used a decimal
approximation (up to 1000 digits) to generate this tree.

Note that this does not have a clearly obvious rhythm, with the childcounts
not appearing to be periodic.

The remaining functions are defined based on a particular $k$-tree, so are
technically also functions of $k$, but for brevity we will often not write the
dependence on $k$ explicitly.

######  Definition 5.

For any $n\in\mathbb{N}$, let
$\text{children}(n)=\\{c\in\mathbb{N}:\lfloor\frac{c}{k}\rfloor=n\\}$ be the
set of children of $n$ in the $k$-tree111Note that technically
$0\in\text{children}(0)$. This makes $h(0)=\lceil k\rceil$, which is
consistent with Lemma 1, as well as Definition 3.. Then
$h(n)=|\text{children}(n)|$ gives the _child-count_ of $n$.

######  Remark 1.

For any $n\in\mathbb{N}_{+}$, $\min(\text{children}(n))=\lceil n\cdot
k\rceil$.

###### Proof.

This follows immediately from the fact that $\lfloor\frac{\lceil
nk\rceil}{k}\rfloor=n$, while $\lfloor\frac{\lceil nk\rceil-1}{k}\rfloor<n$. ∎

We can now formally prove the claim that the average child-count is $k$:

######  Remark 2.

$\lim_{N\to\infty}\frac{1}{N}\sum_{n=0}^{N-1}h(n)=k$.

###### Proof.

Notice that the sum $\sum_{n=0}^{N-1}h(n)$ counts every node from $0$ to the
largest child of the node $N-1$. Thus,

$\sum_{n=0}^{N-1}h(n)=1+\max(\text{children}(N-1))=\min(\text{children}(N))=\lceil
N\cdot k\rceil,$

from Remark 1. Then $\lim_{N\to\infty}\frac{1}{N}\lceil N\cdot k\rceil=k$. ∎

The child count $h(n)\in\\{\lfloor k\rfloor,\lceil k\rceil\\}$, and turns out
to depend only on the quantity $\\{n\cdot k\\}$:

######  Definition 6.

For $n\in\mathbb{N}_{+}$, we call the fractional part $\\{n\cdot k\\}\in[0,1)$
the _count indicator_. The interval $(0,1-\\{k\\}]\subset[0,1)$ is called the
_floor-range_ and its complement $\\{0\\}\cup(1-\\{k\\},1)$ is called the
_ceil-range_.

These definitions are motivated by the following foundational lemma:

######  Lemma 1.

For all $n\in\mathbb{N}$,

$h(n)=\begin{cases}\lfloor k\rfloor&\emph{\text{if }}\\{n\cdot
k\\}\in(0,1-\\{k\\}]\\\ \lceil k\rceil&\emph{\text{otherwise}}.\end{cases}$

###### Proof.

This statement is trivial if $k\in\mathbb{Z}$, so we assume now
$k\notin\mathbb{Z}$ and $\lfloor k\rfloor+1=\lceil k\rceil$.

Notice that $\lfloor\frac{c}{k}\rfloor=n$ iff $c\in[nk,(n+1)k)$, thus $h(n)$
counts the integer points in this interval, whose width is $k=\lfloor
k\rfloor+\\{k\\}$. Intuitively, the first part of the interval $[nk,nk+\lfloor
k\rfloor)$ of width $\lfloor k\rfloor$ will always contain $\lfloor k\rfloor$
integer points. Then when $\\{nk\\}>1-\\{k\\}$, the remaining interval
$[nk+\lfloor k\rfloor,nk+k)$ will “wrap around” one additional integer point.
On the other hand, when $\\{nk\\}<1-\\{k\\}$, the fractional part strictly
increases and does not wrap around an extra integer point.

In the boundary case $\\{nk\\}=0$, $nk\in\mathbb{Z}$, and the interval
contains $\lceil k\rceil$ integer points $nk,nk+1,\ldots,nk+\lfloor k\rfloor$.

In the other boundary case $\\{nk\\}=1-\\{k\\}$, we have $nk+k\in\mathbb{Z}$,
but the interval does not contain this rightmost boundary, and we thus have
$\lfloor k\rfloor$ integer points $nk+k-1,nk+k-2,\ldots,nk+k-\lfloor k\rfloor$
in the interval. ∎

Notice that Lemma 1 implies that for rational $k\in\mathbb{Q}$, the child-
count function $h(n)$ is periodic. See for example Fig. 2. Due to this
periodicity we see that it is a rhythmic tree Definition 3.

Finally, we formally define the row-length sequence:

######  Definition 7.

$(r_{d})_{d=0}^{\infty}$ is the _row-length sequence_ , where $r_{d}$ gives
the number of nodes at depth $d$ in a $k$-tree. More formally, for a node
$n\in\mathbb{N}$, the iterated function sequence $g_{0}=n$ and
$g_{i+1}=\lfloor\frac{g_{i}}{k}\rfloor$ gives the path to the root. Then
$\text{depth}(n)=\max(i:g_{i}=0)$ and
$r_{d}=|\\{n\in\mathbb{N}:\text{depth}(n)=d\\}|$.

The sequence $r_{d}$ then intuitively grows at an exponential rate of $k$ as
how $r_{d}\approx r_{d-1}\cdot k.$ However, due to the rounding off done by
the floor function, it is not exactly this. So, we can define an asymptotics
function on the row lengths.

######  Definition 8.

$\rho(k)=\lim_{d\to\infty}\frac{r_{d}}{k^{d}}$.

Note that we still need to show this limit exists.

##  3\. Asymptotics of $r_{d}$

######  Theorem 1.

For any $k\in(1,\infty)\subset\mathbb{R}$, the constant
$\rho(k)=\lim_{d\to\infty}\frac{r_{d}}{k^{d}}$ exists. Thus
$r_{d}\sim\rho(k)\cdot k^{d}$.

###### Proof.

This is essentially a corollary of Proposition 1 from [3]. To be self-
contained, we produce the proof in its entirety.

Observe that for any $n\in\mathbb{N}_{+}$, $\min(\text{children}(n))=\lceil
n\cdot k\rceil$. We will then consider the sequence $f_{0}=1$ and
$f_{i+1}=\lceil f_{i}\cdot k\rceil$. Notice that, subject to a change in
indexing, this gives the leftmost elements in each row. For example, see

$f_{i}$ | Value
---|---
$f_{0}$ | 1
$f_{1}$ | 2
$f_{2}$ | 4
$f_{3}$ | 7
$f_{4}$ | 12
$f_{5}$ | 20
Table 1. $f_{i}$ values for the $\phi$-tree, shown in Fig. 3.

We thus have $r_{d}=f_{d}-f_{d-1}$ for all $d>0$. Proposition 1 [3] shows
there exists a constant222[3] uses $\alpha$ as the parameter instead of $k$.
$c(k)$ such that $f_{i}\sim c(k)\cdot k^{i}$. This will then imply

$r_{d}\sim c(k)\cdot k^{d}-c(k)\cdot k^{d-1}=\frac{k-1}{k}c(k)\cdot k^{d},$

thus we have $\rho(k)=\frac{k-1}{k}\cdot c(k)$.

To prove Proposition 1 from [3], or that $c(k)$ exists, we consider the
sequence $(\frac{f_{i}}{k^{i}})_{i=0}^{\infty}$. First we show the sequence is
nondecreasing, since

$\frac{f_{i+1}}{k^{i+1}}=\frac{\lceil f_{i}\cdot
k\rceil}{k^{i+1}}\geq\frac{f_{i}}{k^{i}}.$

The sequence is also bounded above. We start with

$\frac{f_{i+1}}{k^{i+1}}=\frac{\lceil f_{i}\cdot
k\rceil}{k^{i+1}}<\frac{f_{i}\cdot
k+1}{k^{i+1}}=\frac{f_{i}}{k^{i}}+\frac{1}{k^{i+1}},$

and then with the base case $\frac{f_{0}}{k^{0}}=1$ we conclude

$\frac{f_{n}}{k_{n}}<\sum_{i=0}^{n}\frac{1}{k^{i}}<\sum_{i=0}^{\infty}\frac{1}{k^{i}}=\frac{1}{1-1/k}=\frac{k}{k-1}.$

Since the sequence $(\frac{f_{i}}{k^{i}})_{i=0}^{\infty}$ is nondecreasing and
bounded from above, the limit $c(k)=\lim_{i\to\infty}(\frac{f_{i}}{k^{i}})$
exists. Moreover, we have the following bounds

$1\leq c(k)\leq\frac{k}{k-1}.$

This proves Proposition 1 from [3], in turn proving the existence of $\rho(k)$
for the asymptotics of $r_{d}$. ∎

######  Corollary 1.

$\frac{k-1}{k}\leq\rho(k)\leq 1$.

###### Proof.

The bounds on $c(k)$ as described in the proof of Theorem 1 alongside the
relationship $\rho(k)=\frac{k-1}{k}\cdot c(k)$ gives the following bounds on
$\rho$. ∎

So, we can now look at approximations of this data. The code for generating
these approximations are attached in Appendix B. Some of the raw data is
contained within Appendix C.

Within [3], Prop 3. they prove that there are jump discontinuities at the
“Josephus points,” or rational numbers of the form $\frac{q}{q-1}$ where
$q\in\mathbb{Z}_{+}.$ This is shown within Fig. 4. Since this is not relevant
to the focus of the paper the proof of this will not be reproduced in its
entirety.

Figure 4. Jump Discontinuities from [3] — In Prop 3. from [3], there are jump
discontinuities of an extremely nice form, where $c(a+\epsilon)=ac(a).$

Zooming out with Fig. 5, we can see what appears to be a more global pattern.

Figure 5. $\rho(k)$ for $1<k\leq 9$ — Approximated $\rho(k)$ for $10^{4}$
points.

One of these patterns appears to be periodic splits within the function. This
is illustrated within Fig. 6. Within $n$ and $(n+1)$ for some integer $n\geq
1$, there are $(n+1)$ periodic ‘visible splits.’ Upon closer inspection, it
appears there are also $(n+1)^{2}$ periodic ‘visible splits,’ with the best
ones happening closest to $(n+1).$ This pattern visually holds, and the splits
get smaller. Upon numerical testing, it is not clear if they share the same
form as those described in [3].

Figure 6. $\rho(k)$ for $1<k\leq 4$ with apparent splits — Approximated
$\rho(k)$ for $10^{4}$ points, and also plotted the points of possible visible
splits that seem to be self similar between $n$ and $n+1$.

##  4\. $k$-Trees with Closed Form $\rho(k)$

Let us then define something that makes it more clear how to describe
indicators.

######  Definition 9.

The child-count indicator graph for some $k>1$ is the set of functions
$\\{f_{i}|i=1,2,\ldots,\lceil k\rceil\\}$, where $f_{i}(x)$ is defined to be
the fractional part of the $i$th child of some $n$ where $x=\\{nk\\},$ so that
$f_{i}:[0,1)\longrightarrow[0,1)$. If $\lceil k\rceil\neq\lfloor k\rfloor$,
then we define $f_{\lceil k\rceil}:(1-\\{k\\},1]\longrightarrow[0,1)$ for the
$ceil{k}$th child of a node $n$ if it exists.

Also, the ceil-range and floor-range are used to say $(1-\\{k\\},1]$ and
$(0,1-\\{k\\}]$ respectively and are helpful to describe the child-count
indicator graph in a more general sense.

What Definition 9 allows us to do is visualize the behavior of the “grand-
children” of a node $n$. If we know it has fractional part $\\{nk\\},$ then we
can see the fractional parts of its children, showing us its children’s
children, leading to the terminology of “grandchildren.”

######  Lemma 2.

Let $k$ satisfy the equation $k=a+\frac{b}{k}$ for $a,b\in\mathbb{Z}$. Then
for any node $n\in\mathbb{N}_{+}$, let $x=\\{n\cdot k\\}\in[0,1)$ be the
count-indicator for $n$. Let $c_{1},\ldots,c_{h(n)}$ be the children of the
node $n$. Then the $i$th smallest child $c_{i}$ has count-indicator

$\\{c_{i}\cdot k\\}=\\{(i-\\{nk\\})\cdot\frac{b}{k}\\}.$

###### Proof.

From Remark 1, we have $c_{1}=\lceil n\cdot k\rceil$, and more generally
$c_{i}=\lceil n\cdot k\rceil+i-1$. We then have count indicator

$\displaystyle\\{c_{i}\cdot k\\}$ $\displaystyle=\\{(\lceil n\cdot
k\rceil+i-1)(a+\frac{b}{k})\\}$ $\displaystyle=\\{(\lceil n\cdot
k\rceil+i-1)\cdot\frac{b}{k}\\}$ $\displaystyle=\\{(n\cdot k+1-\\{n\cdot
k\\}+i-1)\cdot\frac{b}{k})\\}$
$\displaystyle=\\{(i-\\{nk\\})\cdot\frac{b}{k}\\}.$

∎

Some examples of Definition 9 and Lemma 2 are shown in Fig. 10, Fig. 11, Fig.
12, Fig. 13, Fig. 14, located in the Appendix D.

######  Theorem 2.

(Grandparent theorem) Let $a,b\in Z$ such that $a\geq 1$ and $1-a\leq b\leq
a-1$. When $b\geq 0$, then there are always $b$ distinct lines in the ceil-
range of the child count indicator graph for $k=\frac{a+\sqrt{a^{2}+4b}}{2}$.
When $b$ is negative, then there are always $|b|$ distinct lines in the floor-
range of the graph of $k$.

###### Proof.

Let $a,b\in Z$ such that $a\geq 1$ and $1-a\leq b\leq a-1$. Then, let
$k=\frac{a+\sqrt{a^{2}+4b}}{2}$. Using Lemma 2,

$\displaystyle\\{c_{i+1}\cdot k\\}-\\{c_{i}\cdot k\\}$
$\displaystyle=\\{(i+1-x)\cdot\frac{b}{k}\\}-\\{(i-x)\cdot\frac{b}{k}\\}$
$\displaystyle=\\{\frac{b}{k}\\}$
$\displaystyle=\\{k\\},1-\\{k\\}(\text{depending on }b).$

So, the difference between two consecutive childcount indicators is $\\{k\\}$
or $1-\\{k\\}$. The output ceil-range from Lemma 1 is of size $\\{k\\},$ and
the floor-range is then $1-\\{k\\}.$ So, when one childcount indicator falls
above or under that line, the next indicator “flips,” meaning that the amount
of lines in the output ceil-range and floor-range stays constant throughout
$x\in[0,1).$ So, this simplifies the proof in allowing us to choose for $x$ to
show that there are specifically $b$ indicators in the specific output range.

The last step is now to find the number of solutions for $i$ in this
expression, $1-\\{k\\}<\\{\\{k\\}(i-x)\\}<1,$ but now we can choose $x$. If
$x=0,$ then the set of solutions is
$\\{\lfloor\frac{a}{b}\rfloor,\lfloor\frac{2a}{b}\rfloor,\lfloor\frac{3a}{b}\rfloor,\ldots,\lfloor\frac{(b-1)a}{b}\rfloor,\lfloor
a\rfloor\\}$, which contains $b$ solutions (in the nonnegative $b$ case). This
is the set of solutions to $0<\\{\\{k\\}(i-x)\\}<1-\\{k\\}$ for the negative
$b$ case. This then proves the theorem. ∎

######  Theorem 3.

Let $a,b\in\mathbb{Z}$ with $a\geq 1$ and $1-a<b<1+a$. For
$k=\frac{a+\sqrt{a^{2}+4b}}{2}$, the row-length sequence for the $k$-tree
satisfies the linear recurrence

$r_{d}=a\cdot r_{d-1}+b\cdot r_{d-2},$

with base case $r_{0}=1$, $r_{1}=\lceil k\rceil-1$.

###### Proof.

Let $a,b\in\mathbb{Z}$ with $a\geq 1$ and $1-a<b<1+a$. When $b\geq 0,$ we can
see that $\lfloor k\rfloor=a.$ So, $r_{d}\geq ar_{d-1}$, as from Lemma 1 we
know that each element from the previous row contributes at least $\lfloor
k\rfloor$ children. With Theorem 2, we can see that each element in $r_{d-2}$
has $b$ children that have $\lceil k\rceil$ children, which creates the
equality $r_{d}=ar_{d-1}+br_{d-2}.$ Instead, when $b<0,$ then $\lceil
k\rceil=a,$ so $r_{d}\leq ar_{d-1}$. What happens then is $br_{d-2}$ actually
takes away to compensate for the values with $\lfloor k\rfloor$ children. So,

$r_{d}=ar_{d-1}+br_{d-2}.$

∎

######  Corollary 2.

Let $a,b\in\mathbb{Z}$ with $a\geq 1$ and $1-a<b<1+a$. For
$k=\frac{a+\sqrt{a^{2}+4b}}{2}$,

$\rho(k)=\begin{cases}\dfrac{k}{\sqrt{a^{2}+4b}}&b>0\\\\[20.00003pt]
\dfrac{k-1}{\sqrt{a^{2}+4b}}&b<0\\\\[20.00003pt]
\dfrac{a-1}{a}&b=0\end{cases}$

###### Proof.

Because of Theorem 3, we can find a closed formula for $r_{d}$.

As the theorem states, let $a,b\in\mathbb{Z}$ such that $a\geq 1,1-a<b<1+a$,
and let $k_{1}=\frac{a+\sqrt{a^{4}+4b}}{2}$ and
$k_{2}=\frac{a-\sqrt{a^{4}+4b}}{2}$.

Then, for the $k_{1}$ tree,

$r_{d}=ar_{d-1}+br_{d-2}.$

Like any recurrence problem, let $r_{d}=c^{d}$ for some constant $c$. Then,

$\displaystyle c^{d}$ $\displaystyle=ac^{d-1}+bc^{d-2}$ $\displaystyle c^{2}$
$\displaystyle=ac+b$ $\displaystyle c^{2}-ac-b$ $\displaystyle=0$
$\displaystyle c$ $\displaystyle=k_{1},k_{2}.$

So, we can plug in $c$ to $r_{d}$.

$r_{d}=Ak_{1}^{d}+Bk_{2}^{d},$

for some constants $A$ and $B$.

Now, we have to use the starting conditions $r_{0}=1,r_{1}=\lfloor k\rfloor$
to find $A$ and $B$. This is simply done by plugging in $d=0,1$. However, note
that when $b\geq 0$, $\lfloor k\rfloor=a$, and when $b<0$, then $\lfloor
k\rfloor=a-1.$

So, when $b\geq 0$, we can solve two linear systems.

$\begin{bmatrix}1&1&1\\\
k_{1}&k_{2}&a\end{bmatrix}\to\begin{bmatrix}1&0&\dfrac{a+\sqrt{a^{2}+4b}}{2\sqrt{a^{2}+4b}}=\dfrac{k_{1}}{\sqrt{a^{2}+4b}}\\\\[20.00003pt]
0&1&\dfrac{-a+\sqrt{a^{2}+4b}}{2\sqrt{a^{2}+4b}}=\dfrac{-k_{2}}{\sqrt{a^{2}+4b}}\end{bmatrix}$

So, if $b\geq 0$,

$r_{d}=\frac{k_{1}^{d+1}-k_{2}^{d+1}}{\sqrt{a^{2}+4b}}.$

By solving the same original matrix with $a-1$ swapped with $a$, we can get
that when $b<0$,

$r_{d}=\frac{(k_{1}-1)^{d+1}-(k_{2}-1)^{d+1}}{\sqrt{a^{2}+4b}}.$

Now, we can plug these closed formulas into $\rho(k)$. When $b\geq 0,$

$\displaystyle\rho(k)$
$\displaystyle=\lim_{d\to\infty}\dfrac{r_{d}}{k_{1}^{d}}$
$\displaystyle=\lim_{d\to\infty}\dfrac{k_{1}^{d+1}-k_{2}^{d+1}}{k_{1}^{d}\sqrt{a^{2}+4b}}$
$\displaystyle=\lim_{d\to\infty}\dfrac{k_{1}^{d+1}}{k_{1}^{d}\sqrt{a^{2}+4b}}$
$\displaystyle=\dfrac{k_{1}}{\sqrt{a^{2}+4b}}.$

Similarly, when $b<0$,

$\rho(k)=\dfrac{k_{1}-1}{\sqrt{a^{2}+4b}}.$

Note that when $b=0$, then $k=a\in\mathbb{Z}$, so we can also write
$\rho(k)=\dfrac{a-1}{a}$ due to $r_{d}=k^{d-1}(k-1).$ ∎

The values described in Corollary 2 are shown in Table 2.

Table 2. Values of the function $k=\frac{a+\sqrt{a^{2}+4b}}{2}$ for
$a,b\in\mathbb{Z}$. Calculated the formula $k=\dfrac{a+\sqrt{a^{2}+4b}}{2}$
with decimal approximations for $-1\leq a\leq 7,$ and $-6\leq b\leq 8$.

So, by using Corollary 2, we can look at the closed formula for points that we
do know on Fig. 7.

Figure 7. $\rho(k)$ for $1<k\leq 7$ with closed form — Approximated $\rho(k)$
for $10^{4}$ points, and also plotted the points that contain a closed
formula.

##  5\. Conclusion

So, we have found a closed formula for the irrational case, even when the
rational case doesn’t have the same type of behavior. This is a rare case
where irrationality seems to behave more nicely than their rational
counterpart. Also, we have a closed formula for these golden-like trees. This
is interesting and motivates the extension from a $k$-ary tree.

Now, for future work, there seems to be three different modes of progress.
First would simply be to change from $\lfloor\frac{n}{k}\rfloor$ to
$\lceil\frac{n}{k}\rceil$. This would be a $k$-ascending tree.

######  Definition 10.

For any $k\in(1,\infty)\subset\mathbb{R}$, the _$k$ -ascending tree_ is the
rooted tree with nodes in $\mathbb{N}$, where every $n\in\mathbb{N}$ has the
parent $\lceil\frac{n}{k}\rceil$.

Note that the root node of Definition 10 is not going to be $0$.

Figure 8. $\phi$-ascending tree — This generates a $\phi$-ascending tree.
Modified the original code to use the ceiling function instead of floor Figure
9. $\frac{3}{2}$-ascending tree — This generates a $\frac{3}{2}$-ascending
tree. Note that this is still a rhythmic tree with a shift in indexing.

From observation, both examples seem to be shifted (as the root is not $0$),
and Fig. 8 also seems to be flipped. From this it seems that a lot of the same
ideas would carry over very nicely.

Second, it may be interesting to look at more complex $k$’s of the form
$k=a+\frac{b}{k}+\frac{c}{k^{2}}.$ Instead of Theorem 3’s $r_{d}=a\cdot
r_{d-1}+b\cdot r_{d-2}$, we may have the recurrence $r_{d}=a\cdot
r_{d-1}+b\cdot r_{d-2}+c\cdot r_{d-3}.$ It would be nice if this was true, but
it doesn’t seem like it would work due to Theorem 3’s reliance on the
grandparent indicator only going back to $r_{d-2}$.

Also, note that the $\phi$-tree has a regular language that describes its
paths. In [2], they say that all rational numbers $k$ that are not integers do
not have this property in the $k$-tree, so it would be interesting to see
which other irrational numbers have this property. This isn’t true for other
“golden-like” numbers, or numbers of the form $k=a+\frac{b}{k}.$

Furthermore, one of the original reasons that [2] researched rhythms was their
connection to fractional bases. When $k$ is a rational number, the language
describing paths in the $k$-tree form a bijection with the natural numbers in
base $k$. In the irrational $k$ case, the path to a node $n$ in the $k$-tree
could be seen as some sort of base-$k$ representation of $n$, but addition is
not as clear anymore.

Within this paper there was the surprising connection to the Josephus problem
[3], and the article [1] looks specifically at the constant
$\rho(\frac{3}{2})$ in regard to the Collatz Conjecture and mentions the
Josephus problem. The connection to the Collatz Conjecture is light, as one
can simplify the process into this tree. The Collatz Conjecture states that if
you take some $n\in\mathbb{N}$, if it is even divide it by 2, otherwise take
$3n+1$ and continue this process, you will always reach 1. The article talks
about the simplification of this process to just $\lfloor\frac{3n}{2}\rfloor$
to remove the complexity and arrives at the $c(\frac{3}{2})$ constant to
describe the tree that this process generates, which is the same as our
$k$-descending tree.

These connections to the Josephus Problem and others can be further fleshed
out if one found a more direction connection by figuring out how to represent
the Josephus problem with trees. One idea could be that the indices of the
vertices represent the players in the Josephus problem.

This research was presented at the Sofia, Bulgaria WFNMC conference in the
summer of 2022. The authors would like to thank the administration there.

## References

  * [1] James Grime, Kevin Knudson, Pamela Pierce, Ellen Veomett, and Glen Whitney.  Beyond pi and e: a collection of constants.  Math Horizons, 29(1):8–12, 2022.
  * [2] Victor Marsault and Jacques Sakarovitch.  Rhythmic generation of infinite trees and languages.  arXiv preprint arXiv:1403.5190, 2014.
  * [3] Andrew M Odlyzko and Herbert S Wilf.  Functional iteration and the josephus problem.  Glasgow Mathematical Journal, 33(2):235–240, 1991.

Agniv Sarkar Eric Severson

##  Appendix A Tree Generation Code

This is the code used to generate the $k$-tree figures.

⬇

1from decimal import *

2from re import X

3import networkx as nx

4import matplotlib.pyplot as plt

5from math import floor, ceil

6import os

7import random

8

9save_as = 4

10save_location = "./drive/MyDrive/Figures"

11largest_can_draw = 1000

12

13os.makedirs(save_location, exist_ok=True)

14

15def rows(n, k):

16 """ returns r0...r_n in the k-tree """

17 r = row_gen(Decimal(k))

18 return [next(r) for x in range(n)]

19

20def left_gen(k):

21 """ generator for leftmost nodes in the k-tree """

22 yield 0

23 l = 1

24 while True:

25 yield l

26 l = ceil(Decimal(l) * Decimal(k))

27

28def row_gen(k):

29 """ generator for row lengths in the k-tree"""

30 left_iter = left_gen(k)

31 last_left = next(left_iter)

32 while True:

33 next_left = next(left_iter)

34 yield next_left - last_left

35 last_left = next_left

36

37def fancytree(largest, k, save = True):

38 """

39 Draws a k-tree up to largest value given.

40 To use depth instead, pass on the max_n from here:

41

42 depth = 10 # some number of rows

43 max_n = 0

44 row_lengths = rows(depth, k)

45 for row in row_lengths:

46 max_n += row

47

48 """

49 G = nx.DiGraph()

50 fl = []

51 ce = []

52

53 # CONNECTING THE NODES

54 for x in range(1, largest):

55 G.add_edge(floor(x/k), x)

56 if Decimal(k*x) - floor(Decimal(k*x)) < 1 - k + floor(k):

57 fl.append(x)

58 else:

59 ce.append(x)

60

61 # DRAWING THE TREE

62 if (largest > largest_can_draw):

63 return

64 plt.clf()

65 pos = hierarchy_pos(nx.Graph(G), 0)

66 G.add_edge(0, 0)

67

68 # DRAWS THE NODES WITH COLOR

69 nsize = 450

70 lwidth = 3.5

71 nx.draw(G.reverse(), pos, nodelist = [0], node_size = nsize, node_color =
’#ffffff’, node_shape = "o", edgecolors=’#000000’, linewidths = lwidth,
with_labels = True)

72 nx.draw(G.reverse(), pos, label = str(floor(k)) + " Children", nodelist =
fl, node_size = nsize, node_color = ’#ffffff’, node_shape = "s",
edgecolors=’#1357eb’, linewidths = lwidth, with_labels = True)

73 nx.draw(G.reverse(), pos, label = str(ceil(k)) + " Children", nodelist =
ce, node_size = nsize, node_color = ’#ffffff’, node_shape = "o",
edgecolors=’#ab1a1a’, linewidths = lwidth, with_labels = True)

74 if abs(k - int(k)) <= 10**(-5):

75 pass

76 else:

77 plt.legend(loc = ’best’, fontsize = ’xx-large’, markerscale = .45,

78 labelspacing = .7,

79 title_fontsize = ’xx-large’, borderpad = .6, shadow = True)

80 if save:

81 plt.savefig(save_location+str(round(k, save_as))+"tree.pdf", format =
’pdf’, dpi=300)

82 plt.show()

83

84def hierarchy_pos(G, root=None, width=1., vert_gap = 0.2, vert_loc = 0,
leaf_vs_root_factor = 0.5):

85 """

86 Source: https://github.com/ryoji-kubo/CurvatureAnalysis/blob/main/utils.py

87 Source:
https://github.com/BaseMax/BinaryTreeDiagramDrawing/blob/master/tree.py

88 """

89

90 if root is None:

91 if isinstance(G, nx.DiGraph):

92 root = next(iter(nx.topological_sort(G))) #allows back compatibility with
nx version 1.11

93 else:

94 root = random.choice(list(G.nodes))

95

96 if not nx.is_tree(G):

97 raise TypeError(’cannot use hierarchy_pos on a graph that is not a tree’)

98

99 def _hierarchy_pos(G, root, leftmost, width, leafdx = 0.2, vert_gap = 0.2,
vert_loc = 0,

100 xcenter = 0.5, rootpos = None,

101 leafpos = None, parent = None):

102

103 if rootpos is None:

104 rootpos = {root:(xcenter,vert_loc)}

105 else:

106 rootpos[root] = (xcenter, vert_loc)

107 if leafpos is None:

108 leafpos = {}

109 children = list(G.neighbors(root))

110 leaf_count = 0

111 if not isinstance(G, nx.DiGraph) and parent is not None:

112 children.remove(parent)

113 if len(children)!=0:

114 rootdx = width/len(children)

115 nextx = xcenter - width/2 - rootdx/2

116 for child in children:

117 nextx += rootdx

118 rootpos, leafpos, newleaves = _hierarchy_pos(G,child,
leftmost+leaf_count*leafdx,

119 width=rootdx, leafdx=leafdx,

120 vert_gap = vert_gap, vert_loc = vert_loc-vert_gap,

121 xcenter=nextx, rootpos=rootpos, leafpos=leafpos, parent = root)

122 leaf_count += newleaves

123

124 leftmostchild = min((x for x,y in [leafpos[child] for child in children]))

125 rightmostchild = max((x for x,y in [leafpos[child] for child in
children]))

126 leafpos[root] = ((leftmostchild+rightmostchild)/2, vert_loc)

127 else:

128 leaf_count = 1

129 leafpos[root] = (leftmost, vert_loc)

130 return rootpos, leafpos, leaf_count

131

132 xcenter = width/2.

133 if isinstance(G, nx.DiGraph):

134 leafcount = len([node for node in nx.descendants(G, root) if
G.out_degree(node)==0])

135 elif isinstance(G, nx.Graph):

136 leafcount = len([node for node in nx.node_connected_component(G, root) if
G.degree(node)==1 and node != root])

137 rootpos, leafpos, leaf_count = _hierarchy_pos(G, root, 0, width,

138 leafdx=width*1./leafcount,

139 vert_gap=vert_gap,

140 vert_loc = vert_loc,

141 xcenter = xcenter)

142 pos = {}

143 for node in rootpos:

144 pos[node] = (leaf_vs_root_factor*leafpos[node][0] +
(1-leaf_vs_root_factor)*rootpos[node][0], leafpos[node][1])

145 xmax = max(x for x,y in pos.values())

146 for node in pos:

147 pos[node]= (pos[node][0]*width/xmax, pos[node][1])

148 return pos

Listing 1: Tree Generation Code

##  Appendix B $\rho$ Code

This is the code used to generate $\rho(k)$ values.

⬇

1import pandas as pd

2import numpy as np

3from math import floor, ceil, sqrt, log

4from decimal import *

5import matplotlib.pyplot as plt

6

7source = "./drive/MyDrive/Figures/rho.csv"

8

9def resetRho():

10 """ Empties out if error is encountered """

11 df = pd.DataFrame({’val’:[],’out’:[],’precision’:[]})

12 df.to_csv(source, index=False)

13 return df

14

15def getPrec(k, df):

16 """ Returns the stored precision for rho(k) in the csv. 0 if not found """

17 try:

18 return df.loc[min(df.loc[df.val == k].index)].precision

19 except:

20 return 0

21

22def addRho(df, info):

23 """ Adds value to the stored csv. """

24 inside = getPrec(info[0], df) != 0

25 if inside:

26 df.loc[(df["val"] == info[0])] = info

27 else:

28 df.loc[len(df)] = info

29 return df

30

31def writeRho(df):

32 """ Writes to an external csv for storage """

33 df.to_csv(source, index = False)

34

35def left_gen(k):

36 """ Generator for leftmost nodes """

37 yield 0

38 l = 1

39 while True:

40 yield l

41 l = ceil(Decimal(l) * Decimal(k))

42

43def row_gen(k):

44 """ Generator for row lengths """

45 left_iter = left_gen(k)

46 last_left = next(left_iter)

47 while True:

48 next_left = next(left_iter)

49 yield next_left - last_left

50 last_left = next_left

51

52def get_rho(k, min_row_size = 10**3, eps = 10 ** -3, iter_upper_bound = 10
** 2):

53 """

54 If a rho(k) value is not found within rho.csv, generate it.

55 Also runs if a higher precision value

56 precision is defined with:

57 minimum row size to start at, upper bound to stop at, eps, decimal
precision

58 """

59 if k <= 1:

60 return None

61 elif abs(k - int(k)) < 10**(-1000000):

62 return Decimal((k-1))/Decimal(k)

63 row = row_gen(k)

64 last_c = 1

65 r = next(row)

66 c = 1

67 while r < min_row_size:

68 r = next(row)

69 c += 1

70 for i in range(c, iter_upper_bound+c):

71 next_c = Decimal(next(row)) / Decimal(k) ** Decimal(i)

72 if abs(next_c - last_c) < eps:

73 return (next_c + last_c) / 2

74 last_c = next_c

75 return next_c

76

77def find_rho(k, min_row_size = 10**3, eps = 10 ** -3, iter_upper_bound = 10
** 2):

78 """

79 Gets the value of rho(k) to some precision, either through rho.csv or
generating it.

80 """

81 try:

82 df = pd.read_csv(source)

83 except FileNotFoundError:

84 return get_rho(k, min_row_size=min_row_size, eps=eps,
iter_upper_bound=iter_upper_bound)

85 except:

86 df = resetRho()

87 if abs(k - int(k)) < 10**(-100):

88 prec = np.inf

89 else:

90 prec = log(min_row_size*eps**(-1)*iter_upper_bound*getcontext().prec, 10)

91 curPrec = getPrec(k, df)

92 if curPrec < prec or curPrec == 0:

93 rho = get_rho(k, min_row_size=min_row_size, eps=eps,
iter_upper_bound=iter_upper_bound)

94 writeRho(addRho(df, [k, rho, prec]))

95 return rho

96 elif curPrec >= prec:

97 return df.loc[min(df.loc[df.val == k].index)].out

98

99x = np.linspace(1 + 10**(-6), 6, 10**5)

100plt.scatter(x, np.array(list(map(find_rho, x))), s = 0.1)

101plt.xlabel(’k’)

102plt.ylabel(’Ï(k)’)

103plt.title(f’Estimated values of Ï(k) for 100000 random values 1 < k < 6’)

104plt.show()

Listing 2: $\rho(k)$ Approx. Code

##  Appendix C $\rho$ Raw Data

This is some recorded raw data generated by the code. Note that there appears
to be some rounding off that happens as the data is written to the csv file.
The data is actually too large to be typically viewed with tools like
csvsimple. The link is at https://github.com/agniv-the-marker/rho-stuff.

##  Appendix D Grandparent Indicator Figures

Figure 10. $\phi$ Grandparent Indicators Figure 11. $a=3,b=-1$ Grandparent
Indicators Figure 12. $a=3,b=7$ Grandparent Indicators Figure 13. $a=5,b=-3$
Grandparent Indicators Figure 14. $a=5,b=3$ Grandparent Indicators