$\nu^{W}_{3}$ | -2 | -1.98 | 0.09 | 0.08 | -2.42 | -1.6 | -1.99 | 0.06 | 0.06 | -2.25 | -1.72
$\xi(2,1)$ | 0.2 | 0.19 | 0.05 | 0.05 | 0.11 | 0.29 | 0.2 | 0.04 | 0.04 | 0.13 | 0.28
$\xi(3,1)$ | 0.05 | 0.05 | 0.02 | 0.02 | 0.02 | 0.1 | 0.05 | 0.02 | 0.01 | 0.01 | 0.08
$\xi(3,2)$ | 0.5 | 0.49 | 0.07 | 0.06 | 0.4 | 0.6 | 0.5 | 0.05 | 0.04 | 0.41 | 0.6
$\xi(2,3)$ | 0.5 | 0.48 | 0.14 | 0.15 | 0.4 | 0.59 | 0.49 | 0.09 | 0.1 | 0.4 | 0.6
$\beta^{V}_{1}$ | 1 | 0.99 | 0.16 | 0.15 | 0.65 | 1.35 | 0.99 | 0.09 | 0.1 | 0.78 | 1.21
$\beta^{V}_{2}$ | -1 | -1 | 0.13 | 0.14 | -1.2 | -0.83 | -0.99 | 0.09 | 0.1 | -1.11 | -0.88
$\kappa^{V}_{1}$ | 1 | 1.01 | 0.18 | 0.17 | 0.65 | 1.4 | 1 | 0.13 | 0.12 | 0.74 | 1.27
$\kappa^{V}_{2}$ | -1 | -1.01 | 0.29 | 0.28 | -1.64 | -0.44 | -1.01 | 0.21 | 0.21 | -1.52 | -0.57
$\nu^{V}_{1}$ | 1 | 1 | 0.05 | 0.04 | 0.83 | 1.17 | 1 | 0.04 | 0.04 | 0.89 | 1.13
$\nu^{V}_{2}$ | 1 | 1.01 | 0.06 | 0.05 | 0.82 | 1.21 | 1 | 0.04 | 0.03 | 0.87 | 1.14
$\nu^{V}_{3}$ | -2 | -1.99 | 0.07 | 0.05 | -2.29 | -1.73 | -2 | 0.04 | 0.04 | -2.2 | -1.81
Table 3: Summary statistics of the parameter estimates for the Mreps
replications in the simulation runs for sample sizes $n=\\{50,100\\}$. The
columns are the true values of the model parameters, the sample mean of the
estimates, the sample standard deviations, the asymptotic standard errors, the
2.5% percentile, and the 97.5% percentile. The sample standard deviations are
estimates of the standard errors of the estimators. The asymptotic standard
errors are the means of the asymptotic standard errors over the Mreps
replications. The RCR component includes model parameters $\Lambda$s,
$\alpha$s and those parameters with superscript $R$. The LM component includes
model parameters $\eta$s and those parameters with superscript $W$. The HS
component includes model parameters $\xi$s and those parameters with
superscript $V$.

The graphical summary of these centered estimates for Mreps replications of
$n=50$ units is given in Figure 5. Centering for each estimator is done by
subtracting the true value of the parameter being estimated. We observe that
the medians of all these centered estimates are close to $0$. We also observe
some outliers, but most of the centered Mreps estimates are close to $0$. In
Figure 6, we show three types of plots for the baseline survivor function for
each risk in the RCR component. The true Weibull type baseline survivor
function is plotted in red color, the overlaid plots of a random selection of
ten estimates of the baseline survivor functions are in green color, and the
mean baseline survivor function based on the Mreps estimates is shown in blue
color. We observe that there is close agreement between the true (red) and
mean (blue) curves. Based on these simulation studies, the semi-parametric
estimation procedure for the joint model appears to provide reasonable
estimates of the true finite- and infinite-dimensional model parameters, at
least for the choices of parameter values for these particular simulations.
Further simulation and analytical studies are still needed to substantively
assess the performance of the semi-parametric estimation procedure for the
proposed joint model.

---
Figure 5: Boxplots of the centered parameter estimates from Mreps replications
for simulated data sets each with $n=50$ units. Centering is done by
subtracting the true parameter value from each of the Mreps estimates.
---
Figure 6: Overlaid plots of the true baseline survivor function (in red), ten
simulated estimates of the baseline survivor function (in green), and the mean
baseline survivor function based on Mreps simulations (in blue) for each of
the three risks in the RCR component.

##  8 Illustration on a Real Data Set

To illustrate our estimation procedure on a real data set, we apply the joint
model and the semi-parametric estimation procedure to a medical data set with
$n=150$ patients diagnosed with metastatic colorectal cancer which cannot be
controlled by curative surgeries. This data set was gathered in France from
2002–2007 and was used in [10]. It consists of two data sets which are
deposited in the frailtypack package in the R Library [26]: data set
colorectal.Longi and data set colorectal. The data set colorectal.Longi
includes the follow-up period, in years, of the patient’s tumor size
measurements. The times of first measurements of tumor size vary from patient
to patient, so to have all of them start at time ‘zero’, our artificial time
origin, we consider these first measurements as their initial states.
Subsequent times of measuring tumor size are then in terms of the lengths of
time from their time origin. There were a total of $906$ tumor size
measurements for all the patients. In order to conform to our discrete-valued
LM model, we classify (arbitrarily) the tumor size into three categories
(states): $1$, $2$, and $3$, if the tumor size belongs in the intervals
$[3.4,6.6]$, $[2,3.4]$, and $(0,2)$, respectively. Since the tumor size is
only measured at discrete times, instead of continuously, we assumed that the
tumor size state is constant between tumor size measurement times, and
consequently the tumor size process could only transition at the times in
which tumor size is measured. This assumption is most likely unrealistic, but
we do so for the purpose of illustration. The data set colorectal contains
some information about the patient’s ID number, covariates $X_{1}$ and
$X_{2}$, with $X_{1}=1(0)$ if patient received treatment C (S); $X_{2}$
consists of two dummy variables, with $X_{2}=(0,0),(1,0),(0,1)$ if the initial
WHO performance status is $0,1,2$, respectively; the time (in years) of each
occurrence of a new lesion since baseline measurement time; and the final
right-censored or death time. There were $289$ occurrences of new lesions and
$121$ patients died during the study.

Clearly, this data set is a special case of the type of data appropriate for
our proposed joint model, having only one type of recurrent event, one
absorbing health status state (dead), and one transient health status state
(alive). We assume the effective age
$\mathcal{E}_{i}(s)=s-S_{iN_{i}^{R}(s-)}^{R}$ after each occurrence of a new
lesion, and use $\rho(k);\alpha)=\alpha^{\log(1+k)}$. The unknown model
parameters in the RCR (here, just a recurrent event) component of the model
includes $\alpha$ in the $\rho(\cdot)$ function,
$\beta^{R}=[\beta^{R}_{1},\beta^{R}_{2},\beta^{R}_{3}]$ for the covariates,
and $[\kappa^{R}_{1},\kappa^{R}_{2}]$ for the LM state. The unknown model
parameters in the LM component of the model are $\beta^{W}$ for the covariates
and $\nu^{W}_{1}$ for the recurrent event process. The unknown parameters in
the HS component of the model includes $\beta^{V}$ for the covariates,
$[\kappa^{V}_{1},\kappa^{V}_{2}]$ for the LM state, and $\nu^{V}_{1}$ for the
recurrent event counting process.

We fitted the joint model to this data set. The resulting model parameter
estimates along with the information-based standard error estimates are given
in Table 4. The standard errors are obtained by taking square roots of the
diagonal elements of the observed inverse of the profile likelihood
information matrix. Based on these estimates, we could also perform hypothesis
tests. Thus, for instance, the $p$-values associated with the two-tailed
hypothesis tests are also given in the table. The null hypothesis being tested
for $\alpha$ is that $H_{0}:\alpha=1$, while the null hypotheses for the other
model parameters are that their true parameter values are zeros. We test
$\alpha=1$ instead of $\alpha=0$ because $\alpha=1$ means that the
accumulating number of recurrent event occurrences does not have an impact in
subsequent recurrent event occurrences. From the values in Table 4, the
estimate of $\alpha$ is less than one, which may indicate that each occurrence
of new lesion decreases the risk of future occurrences of new lesions, though
from the result of the statistical test we cannot conclude statistically that
$\alpha<1$. Based on the set of $p$-values, we find that the initial WHO
performance state of $1$ and the tumor size state of $2$ are associated with
decreased risk of new lesion occurrences, while an initial WHO performance
state of $2$ is associated with an increased risk of new lesion occurrences.
An initial WHO performance state of $2$ and the number of occurrences of new
lesions are associated with an increased risk of death in the health status.

Finally, we want to emphasize the importance of the effective age process. An
inappropriate effective age may lead to misleading estimates. Parameter model
estimates under a mis-specified effective age could lead to biases and
potentially misleading conclusions. This is one aspect where domain
specialists and statisticians need to consider when assessing the impact of
interventions since the specification of the effective or virtual age have
important and consequential implications. For more discussions about effective
or virtual ages, see the recent papers [12, 4], the last one also touching on
the situation where the virtual age function depends on unknown parameters.

Parameter | Estimate | Est. Standard Error |  $p$-value
---|---|---|---
$\alpha$ | 0.77 | 0.22 | 0.30
$\beta^{R}_{1}$ | -0.16 | 0.20 | 0.43
$\beta^{R}_{2}$ | -0.42 | 0.21 | 0.05
$\beta^{R}_{3}$ | 0.88 | 0.41 | 0.03
$\kappa^{R}_{1}$ | -0.52 | 0.27 | 0.05
$\kappa^{R}_{2}$ | -0.25 | 0.25 | 0.31
$\beta^{W}_{1}$ | 0.08 | 0.25 | 0.75
$\beta^{W}_{2}$ | -0.00 | 0.25 | 0.99
$\beta^{W}_{3}$ | -0.51 | 0.73 | 0.48
$\nu^{W}_{1}$ | -0.23 | 0.18 | 0.19
$\beta^{V}_{1}$ | 0.49 | 0.30 | 0.10
$\beta^{V}_{2}$ | -0.11 | 0.30 | 0.71
$\beta^{V}_{3}$ | 1.20 | 0.41 | 0.00
$\kappa^{V}_{1}$ | -0.43 | 0.31 | 0.16
$\kappa^{V}_{2}$ | -0.18 | 0.33 | 0.59
$\nu^{V}_{1}$ | 0.71 | 0.14 | 0.00
Table 4: Parameter estimates, information-based standard errors, and
$p$-values for the RCR, LM, and HS processes based on the real data set. The
$p$-value is based on the two-tailed hypothesis test that the model parameter
is zero (except for the parameter $\alpha$ where $H_{0}:\alpha=1$). The top
block includes the model parameters in the RCR component, the middle block
includes the model parameters in the LM component, and the bottom block
includes the model parameters in the HS component.

##  9 Concluding Remarks

For the general class of joint models for recurrent competing risks,
longitudinal marker, and health status proposed in this paper, which
encompasses many existing models considered previously, there are still
numerous aspects that need to be addressed in future studies. Foremost among
these aspects is a more refined analytical study of the finite-sample and
asymptotic properties of the estimators of model parameters, together with
other inferential and prediction procedures. The finite-sample and asymptotic
results could be exploited to enable performing tests of hypothesis and
construction of confidence regions for model parameters. There is also the
interesting aspect of computationally estimating the standard errors of the
estimators. How would a bootstrapping approach be implemented in this
situation? Another important problem that needs to be addressed is how to
perform goodness-of-fit and model validation for this joint model. Though the
class of models is very general, there are still possibilities of model mis-
specifications, such as, for example, in determining the effective age
processes, or in the specification of the $\rho_{q}(\cdot)$-functions. What
are the impacts of such model mis-specifications? Do they lead to serious
biases that could potentially result in misleading conclusions? These are some
of the problems whose solutions await further studies.

A potential promise of this joint class of models is in precision medicine.
Because all three components (RCR, LM, HS) are taken into account
simultaneously, in contrast to a marginal modeling approach, the synergy that
this joint model allows may improve decision-making – for example, in
determining interventions to be performed for individual units. In this
context, it is of utmost importance to be able to predict in the future the
trajectories of the HS process given information at a given point in time
about all three processes. Thus, an important problem to be dealt with in
future work is the problem of forecasting using this joint model. How should
such forecasting be implemented? This further leads to other important
questions. One is determining the relative importance of each of the
components in this prediction problem. Could one ignore other components and
still do as well relative to a joint model-based prediction approach? If there
are many covariates, how should the important covariates among these numerous
covariates be chosen in order to improve prediction of, say, the time-to-
absorption?

Finally, though our class of joint models is a natural extension of earlier
models dealing with either recurrent events, competing recurrent events,
longitudinal marker, and terminal events, one may impugn it as not realistic,
but instead view it as more of a futuristic class of models, since existing
data sets were not gathered in the manner for which these joint models apply.
For instance, in the example pertaining to gout in Section 2, the SUR level
and CKD status are not continuously monitored. However, with the advent of
smart devices, such as smart wrist watches, embedded sensors, black boxes,
etc., made possible by miniaturized technology, high-speed computing, almost
limitless cloud-based memory capacity, and availability of rich cloud-based
databases, the era is, in our opinion, fast approaching when continuous
monitoring of longitudinal markers, health status, occurrences of different
types of recurrent events, be it on a human being, an experimental animal or
plant, a machine such as an airplane or car, an engineering system, a business
entity, etc. will become more of a standard rather than an exception. By
developing the models and methods of analysis for such future complex and
advanced data sets, even before they become available and real, will hasten
and prepare us for their eventual arrival.

## Acknowledgments

Some materials in this work are based on L. Tong’s PhD dissertation research
at the University of South Carolina. E. Peña is currently Program Director in
the Division of Mathematical Sciences (DMS) at the National Science Foundation
(NSF). As a consequence of this position, he receives support for research,
which included work in this manuscript, under NSF Grant 2049691 to the
University of South Carolina. Any opinions, findings, and conclusions or
recommendations expressed in this material are those of the authors and do not
necessarily reflect the views of NSF. P. Liu acknowledges the 2017–2019 Summer
Research Grants from Bentley University. E. Peña acknowledges NSF Grant
2049691 and NIH Grant P30 GM103336-01A1.

## References

  * [1] P. K. Andersen, Ø. Borgan, R. D. Gill, and N. Keiding.  Statistical models based on counting processes.  Springer Series in Statistics. Springer-Verlag, New York, 1993.
  * [2] P. K. Andersen and R. D. Gill.  Cox’s regression model for counting processes: a large sample study.  Ann. Statist., 10(4):1100–1120, 1982.
  * [3] R. E. Barlow and F. Proschan.  Statistical theory of reliability and life testing.  Holt, Rinehart and Winston, Inc., New York-Montreal, Que.-London, 1975\.  Probability models, International Series in Decision Processes, Series in Quantitative Methods for Decision Making.
  * [4] E. Beutner.  A review of effective age models and associated non- and semiparametric methods.  To appear in Econometrics and Statistics, 2021.
  * [5] P. Blanche, C. Proust-Lima, L. Loubère, C. Berr, J.-F. Dartigues, and H. Jacqmin-Gadda.  Quantifying and comparing dynamic predictive accuracy of joint models for longitudinal marker and time-to-event in presence of censoring and competing risks.  Biometrics, 71(1):102–113, 2015.
  * [6] O. Borgan.  Maximum likelihood estimation in parametric counting process models, with applications to censored failure time data.  Scandinavian Journal of Statistics, 11:1–16, 1984.
  * [7] Q. Cai, M.-C. Wang, and K. C. G. Chan.  Joint modeling of longitudinal, recurrent events and failure time data for survivor’s population.  Biometrics, 73(4):1150–1160, 2017.
  * [8] D. R. Cox.  Regression models and life-tables.  J. Roy. Statist. Soc. Ser. B, 34:187–220, 1972.
  * [9] C. Dorado, M. Hollander, and J. Sethuraman.  Nonparametric estimation for a general repair model.  Ann. Statist., 25(3):1140–1160, 1997.
  * [10] M. Ducreux, D. Malka, J. Mendiboure, P.-L. Etienne, P. Texereau, D. Auby, P. Rougier, M. Gasmi, M. Castaing, M. Abbas, P. Michel, D. Gargot, A. Azzedine, C. Lombard-Bohas, P. Geoffroy, B. Denis, J.-P. Pignon, L. Bedenne, O. Bouché, and F. F. de Cancérologie Digestive (FFCD) 2000-05 Collaborative Group.  Sequential versus combination chemotherapy for the treatment of advanced colorectal cancer (ffcd 2000-05): an open-label, randomised, phase 3 trial.  The Lancet. Oncology, 12(11):1032–1044, Oct 2011.
  * [11] A. Efendi, G. Molenberghs, E. N. Njagi, and P. Dendale.  A joint model for longitudinal continuous and time-to-event outcomes with direct marginal interpretation.  Biom. J., 55(4):572–588, 2013.
  * [12] M. Finkelstein and J. H. Cha.  Virtual age, is it real? Discussing virtual age in reliability context.  Applied Stochastic Models in Business and Industry, 37(1):3–16, January 2021.
  * [13] H. K. Gjessing, K. Røysland, E. A. Peña, and O. O. Aalen.  Recurrent events and the exploding Cox model.  Lifetime Data Anal., 16(4):525–546, 2010.
  * [14] J. R. González, E. A. Peña, and E. H. Slate.  Modelling intervention effects after cancer relapses.  Statistics in medicine, 24(24):3959–3975, 2005.
  * [15] J. Han, E. H. Slate, and E. A. Peña.  Parametric latent class joint model for a longitudinal biomarker and recurrent events.  Stat. Med., 26(29):5285–5302, 2007.
  * [16] R. Henderson, P. Diggle, and A. Dobson.  Joint modelling of longitudinal measurements and event time data.  Biostatistics, 1(4):465–480, 2000.
  * [17] S. Kim, D. Zeng, L. Chambless, and Y. Li.  Joint models of longitudinal data and recurrent events with informative terminal event.  Statistics in Biosciences, 4(2):262–281, 2012.
  * [18] A. Krol, A. Mauguen, Y. Mazroui, A. Laurent, S. Michiels, and V. Rondeau.  Tutorial in joint modeling and prediction: a statistical software for correlated longitudinal outcomes, recurrent events and a terminal event.  Journal of Statistical Software, Articles, 81(3):1–52, 2017.
  * [19] L. Liu, R. A. Wolfe, and X. Huang.  Shared frailty models for recurrent events and a terminal event.  Biometrics, 60(3):747–756, 2004.
  * [20] P. Liu and E. Peña.  Dynamic modeling & analysis of recurrent competing risks and a terminal event.  Statistical, Stochastic and Data Analysis Methods & Applications, 1(1):111–127, 2015.
  * [21] A. Mauguen, B. Rachet, S. Mathoulin-Pélissier, G. MacGrogan, A. Laurent, and V. Rondeau.  Dynamic prediction of risk of death using history of cancer recurrences in joint frailty models.  Statistics in Medicine, 32(30):5366–5380, 2013.
  * [22] A. C. McLain, R. Sundaram, and G. M. Buck Louis.  Joint analysis of longitudinal and survival data measured on nested timescales by using shared parameter models: an application to fecundity data.  J. R. Stat. Soc. Ser. C. Appl. Stat., 64(2):339–357, 2015.
  * [23] E. A. Peña.  Dynamic modeling and statistical analysis of event times.  Statist. Sci., 21(4):487–500, 2006.
  * [24] E. A. Peña.  Asymptotics for a class of dynamic recurrent event models.  J. Nonparametr. Stat., 28(4):716–735, 2016.
  * [25] E. A. Peña, E. H. Slate, and J. R. González.  Semiparametric inference for a general class of models for recurrent events.  Journal of Statistical Planning and Inference, 137(6):1727–1747, 2007.
  * [26] R Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2018.
  * [27] X. Song, M. Davidian, and A. A. Tsiatis.  A semiparametric likelihood approach to joint modeling of longitudinal and time-to-event data.  Biometrics, 58(4):742–753, 2002.
  * [28] L. L. Taylor and E. A. Peña.  Nonparametric estimation with recurrent competing risks data.  Lifetime Data Analysis, 20(4):514–537, 2014.
  * [29] A. Tsiatis, V. DeGruttola, and M. Wulfsohn.  Modeling the relationship of survival to longitudinal data measured with error. applications to survival and cd4 counts in patients with aids.  Journal of the American Statistical Association, 90(429):27–37, 1995.
  * [30] A. A. Tsiatis and M. Davidian.  Joint modeling of longitudinal and time-to-event data: an overview.  Statistica Sinica, 14(3):809–834, 2004.
  * [31] T. Velavan, S. Kuk, L. Linh, and et al.  Longitudinal monitoring of laboratory markers characterizes hospitalized and ambulatory covid-19 patients.  Sci Rep, 11 (14471), 2021.
  * [32] M. S. Wulfsohn and A. A. Tsiatis.  A joint model for survival and longitudinal data measured with error.  Biometrics, 53:330–339, 1997.


[108, 111], and the middle term represents the effects of noise-induced
selection. Note that in finite populations, mean fitness can change
directionally even when there is no standing variation in fitness
($\sigma^{2}_{w}=0$) as long as
$\sigma^{2}_{b^{\textrm{(ind)}}}\neq\sigma^{2}_{d^{\textrm{(ind)}}}$. Equation
S61 has been independently derived for discrete traits in a recent preprint
focused on life-history evolution [112].

### The infinite population limit

In the infinite population limit ($K\to\infty$), the expectations in equation
Eq. S61 become superfluous (since the process is now deterministic), and we
thus obtain

$\frac{d\overline{w}}{dt}=\sigma^{2}_{w}+\overline{\hphantom{a}\frac{\partial
w}{\partial t}\hphantom{a}}$ (S62)

which is Fisher’s fundamental theorem in the presence of ecological changes to
fitness [108, 111].

##  S6 Deriving a stochastic version of gradient dynamics

In this section, I derive an equation of gradient dynamics for finite,
fluctuating equations. I begin by once more listing my assumptions. I assume

  * •

Rare mutations, _i.e._ $\mu$ is infinitesimally small.

  * •

Small mutational effects with ‘almost faithful’ reproduction, meaning
$Q(x|\phi)$ is infinitesimally small.

  * •

Strong selection, meaning that variants with low relative fitness are
immediately eliminated from the population and the population remains sharply
peaked around a small number of trait values. Mathematically, this means that
if we begin with a monomorphic population $\phi(x,0)=N_{K}(0)\delta_{y_{0}}$
for some constants $N_{K}(0)>0$ and $y_{0}\in\mathcal{T}$, the population
remains sufficiently well clustered for some time $t>0$ that we can continue
to approximate the distribution $\phi(x,t)$ as a Dirac Delta mass
$N_{K}(t)\delta_{y(t)}$ that is moving across the trait space according to a
trajectory dictated by a function $y(t)$ (to be found). Note that we are not
assuming that the population is exactly given by a delta mass, only that it is
sufficiently sharply peaked that it can be approximated as one (i.e. ‘almost’
every individual has trait value $y$).

I will call these together the strong-selection-weak-mutation (SSWM)
assumption. Let $\int=\int_{\mathcal{T}}$ for notational convenience. Under
SSWM, we wish to obtain an equation for $y(t)$. Note first that we have:

$p(x,t)=\frac{\phi(x,t)}{N_{K}(t)}=\frac{N_{K}(t)\delta_{y(t)}}{N_{K}(t)}=\delta_{y}(t)$
(S63)

Further, for any type level quantity $f(x|\phi)$, we have

$\overline{f}(t)=\int f(x|\phi)p(x,t)dx=\int
f(x|N_{K}\delta_{y(t)})\delta_{y(t)}dx=f(y(t)|N_{K}\delta_{y(t)})$ (S64)

We can now substitute $f(x|\phi)=x$ into our stochastic Price equation (Eq.
S55) and, by virtue of our SSWM assumptions, neglect all terms of the form
$\mu Q(x|\phi)$. Upon doing this and using the expressions for $\overline{x}$
and $p(x,t)$ from above, we obtain

$\frac{dy}{dt}=\textrm{Cov}(w,x)-\frac{1}{KN_{K}(t)}\textrm{Cov}(\tau,x)+\frac{1}{\sqrt{K}N_{K}}\frac{dW_{y}}{dt}$
(S65)

where

$\frac{dW_{y}}{dt}=\int\limits_{\mathcal{T}}(x-y(t))\sqrt{\tau(x|N_{K}(t)\delta_{y}(t))}\dot{W}(x,t)dx$
(S66)

Let $f(x|\phi)$ be any type-level quantity. We can calculate the statistical
covariance of $f$ with $x$ under SSWM as

$\displaystyle\textrm{Cov}(f,x)$
$\displaystyle=\int(f(x|\phi)-\overline{f})(x-\overline{x})p(x,t)dx$ (S67)
$\displaystyle=\int(f(x|N_{K}\delta_{y(t)})-f(y(t)|N_{K}\delta_{y(t)}))(x-y(t))p(x,t)dx$
(S68)

Now, since under SSWM the population at time $t$ is sharply peaked around
$y(t)$, we only need to worry about the value of $f(x|N_{K}\delta_{y(t)})$ at
trait values that are (infinitesimally) close to $y(t)$. We can thus Taylor
expand $f$ as:

$f(x|N_{K}\delta_{y(t)})=f(y(t)|N_{K}\delta_{y(t)})+(x-y)\frac{\partial
f}{\partial z}(z|N_{K}\delta_{y(t)})\bigg{|}_{z=y}+\cdots$ (S69)

and neglect all higher order terms. Using the relation Eq. S69 in Eq. S68, we
obtain

$\displaystyle\textrm{Cov}(f,x)$ $\displaystyle=\int(x-y)^{2}\frac{\partial
f}{\partial z}(z|N_{K}\delta_{y(t)})\bigg{|}_{z=y}p(x,t)dx$ (S70)
$\displaystyle=\int(x-y)^{2}p(x,t)dx\frac{\partial f}{\partial
z}(z|N_{K}\delta_{y(t)})\bigg{|}_{z=y}$ (S71)
$\displaystyle=\int(x-\overline{x})^{2}p(x,t)dx\frac{\partial f}{\partial
z}(z|N_{K}\delta_{y(t)})\bigg{|}_{z=y}$ (S72)
$\displaystyle\Rightarrow\textrm{Cov}(f,x)$
$\displaystyle=\sigma^{2}_{x}(t)\frac{\partial f}{\partial
z}(z|N_{K}\delta_{y(t)})\bigg{|}_{z=y}$ (S73)

Finally, using Eq. S73 in Eq. S65 with $f=w$ and $f=\tau$ yields

$\frac{dy}{dt}=\sigma^{2}_{x}(t)\frac{\partial}{\partial
z}\left[w(z|N_{K}\delta_{y(t)})-\frac{1}{KN_{K}(t)}\tau(z|N_{K}\delta_{y(t)})\right]\bigg{|}_{z=y(t)}+\frac{1}{\sqrt{K}N_{K}}\frac{dW_{y}}{dt}$
(S74)

which is equation 30 in the main text. Note that strictly speaking, if
$\phi(x,t)=N_{K}\delta_{y(t)}$ exactly, then $\sigma^{2}_{x}\equiv 0$. This
just reflects our assumption that mutations are vanishingly rare and mutants
are sampled from infinitesimally close to the resident value, _i.e._ that the
population evolves in very small mutational steps. Thus, for our equation to
make sense, the limits ${\mu\to 0}$, ${Q(x|\phi)\to 0}$, and
${\phi(x,t)\to\delta_{y(t)}N_{K}(t)}$ must be taken simultaneously such that
$\sigma^{2}_{x}(t)\not\to 0$, a standard assumption in such heuristic
derivations of gradient dynamics [105, 113]. More detailed mathematical
arguments are required to ensure that this limit ‘makes sense’. Such
convergence has been proved rigorously for the standard gradient equation
using sophisticated mathematical tools grounded in martingale theory [99].

##  S7 Detecting phenotypic clustering and speciation using spectral methods

In this section, I present a general analytical technique to study phenotypic
clustering using the Fourier transform. To do this in full generality, we
first require an approximation known as the weak noise approximation.

###  S7.1 The weak noise approximation

If stochasticity takes the form of weak fluctuations about a deterministic
trajectory, we can carry out a functional analog of the ‘linear noise
approximation’ or ‘weak noise approximation’ [109] to obtain a linear
functional Fokker-Planck equation for the population density field. In the
next section, we will use this linear approximation to illustrate how
phenotypic clustering can be detected.

Assume that $\psi(x,t)$ is the deterministic trajectory of the density field
in the infinite population limit, obtained as the solution to Eq. 19. Consider
a new process $\\{\zeta(\cdot,s)\\}_{s\geq 0}$ which measures the fluctuations
of $\phi(x,t)$ from the deterministic trajectory $\psi(x,t)$. More precisely,
let us introduce the new variables:

$\displaystyle\zeta(x,s)$ $\displaystyle=\sqrt{K}(\phi(x,t)-\psi(x,t))$ (S75)
$\displaystyle s$ $\displaystyle=t$ $\displaystyle\tilde{P}(\zeta,s)$
$\displaystyle=\frac{1}{\sqrt{K}}P(\phi,t)$

where we have introduced a new time variable $s=t$ just to makes the
substitutions/calculations cleaner. Note that the following relations hold:

$\displaystyle\frac{\delta F[\zeta]}{\delta\phi(x)}$
$\displaystyle=\int\limits_{\mathcal{T}}\frac{\delta
F[\zeta]}{\delta\zeta(y)}\frac{\delta\zeta(y)}{\delta\phi(x)}dy=\sqrt{K}\frac{\delta
F[\zeta]}{\delta\zeta(x)}\ $ (S76) $\displaystyle\frac{\partial}{\partial s}$
$\displaystyle=\frac{\partial}{\partial t}$ (S77)

Furthermore, for any $\zeta\in\mathcal{M}_{K}(\mathcal{T})$, we have:

$\displaystyle\frac{\partial\tilde{P}}{\partial t}(\zeta,s)$
$\displaystyle=\frac{\delta\tilde{P}}{\delta\zeta}\frac{\partial\zeta}{\partial
t}+\frac{\partial\tilde{P}}{\partial s}\frac{\partial s}{\partial t}$
$\displaystyle=\frac{\delta\tilde{P}}{\delta\zeta}\left(-\sqrt{K}\frac{\partial\psi}{\partial
t}\right)+\frac{\partial\tilde{P}}{\partial s}$
$\displaystyle=-\sqrt{K}\frac{\delta}{\delta\zeta}\\{\mathcal{A}^{-}(x|\psi)\tilde{P}(\zeta,s)\\}+\frac{\partial\tilde{P}}{\partial
s}$ (S78)

Reformulating equation Eq. 15 in terms of the new variables Eq. S75 and using
the relations Eq. S76, Eq. S77 and Eq. S78, we obtain:

$\begin{split}-\sqrt{K}\frac{\delta}{\delta\zeta(x)}\\{\mathcal{A}^{-}(x|\psi)\tilde{P}(\zeta,s)\\}+\frac{\partial\tilde{P}}{\partial
s}=\int\limits_{\mathcal{T}}\left[-\left(\sqrt{K}\frac{\delta}{\delta\zeta(x)}\right)\\{\mathcal{A}^{-}\left(x\bigg{|}\psi+\frac{\zeta}{\sqrt{K}}\right)\tilde{P}(\zeta,s)\\}\right]dx\\\
+\int\limits_{\mathcal{T}}\left[\frac{1}{2K}\left(K\frac{\delta^{2}}{\delta\zeta(x)^{2}}\right)\\{\mathcal{A}^{+}\left(x\bigg{|}\psi+\frac{\zeta}{\sqrt{K}}\right)\tilde{P}(\zeta,s)\\}\right]dx\end{split}$

and rearranging gives us:

$\displaystyle\begin{split}\frac{\partial\tilde{P}}{\partial
s}&=-\sqrt{K}\int\limits_{\mathcal{T}}\frac{\delta}{\delta\zeta(x)}\left\\{\left(\mathcal{A}^{-}\left(x\bigg{|}\psi+\frac{\zeta}{\sqrt{K}}\right)-\mathcal{A}^{-}(x|\psi)\right)\tilde{P}(\zeta,s)\right\\}dx\\\
&+\frac{1}{2}\int\limits_{\mathcal{T}}\frac{\delta^{2}}{\delta\zeta(x)^{2}}\\{\mathcal{A}^{+}\left(x\bigg{|}\psi+\frac{\zeta}{\sqrt{K}}\right)\tilde{P}(\zeta,s)\\}dx\end{split}$
(S79)

We will now Taylor expand our functionals about $\psi$ (I assume that this is
possible). Thus, we have the expansions:

$\displaystyle\mathcal{A}^{-}\left(x\bigg{|}\psi+\frac{\zeta}{\sqrt{K}}\right)$
$\displaystyle=\mathcal{A}^{-}\left(x|\psi\right)+\frac{1}{\sqrt{K}}\int\limits_{\mathcal{T}}\zeta(y)\frac{\delta}{\delta\psi(y)}\\{\mathcal{A}^{-}(y|\psi)\\}dy+\cdots$
$\displaystyle\mathcal{A}^{+}\left(x\bigg{|}\psi+\frac{\zeta}{\sqrt{K}}\right)$
$\displaystyle=\mathcal{A}^{+}\left(x|\psi\right)+\frac{1}{\sqrt{K}}\int\limits_{\mathcal{T}}\zeta(y)\frac{\delta}{\delta\psi(y)}\\{\mathcal{A}^{+}(y|\psi)\\}dy+\cdots$

I also assume that $\tilde{P}$ can be expanded as

$\displaystyle\tilde{P}$
$\displaystyle=\sum\limits_{n=0}^{\infty}\tilde{P}_{n}\left(\frac{1}{\sqrt{K}}\right)^{n}$

substituting these expansions into equation Eq. S79, equating coefficients of
powers of $1/K$, and truncating at the lowest order term, we have:

$\displaystyle\frac{\partial\tilde{P}_{0}}{\partial
s}(\zeta,s)=\int\limits_{\mathcal{T}}\left[-\frac{\delta}{\delta\zeta(x)}\left\\{\int\limits_{\mathcal{T}}\zeta(y)\frac{\delta}{\delta\psi(y)}\\{\mathcal{A}^{-}(y|\psi)\\}dy\tilde{P}_{0}(\zeta,s)\right\\}+\frac{1}{2}\mathcal{A}^{+}(x|\psi)\frac{\delta^{2}}{\delta\zeta(x)^{2}}\\{\tilde{P}_{0}(\zeta,s)\\}\right]dx$

We thus arrive at the functional Fokker-Planck equation:

$\frac{\partial\tilde{P}_{0}}{\partial
s}(\zeta,s)=\int\limits_{\mathcal{T}}\left(-\frac{\delta}{\delta\zeta(x)}\left\\{\mathcal{D}_{\zeta}[\mathcal{A}^{-}](x)\tilde{P}_{0}(\zeta,s)\right\\}+\frac{1}{2}\mathcal{A}^{+}(x|\psi)\frac{\delta^{2}}{\delta\zeta(x)^{2}}\\{\tilde{P}_{0}(\zeta,s)\\}\right)dx$
(S80)

where

$\mathcal{D}_{\zeta}[\mathcal{A}^{-}](x)=\int\limits_{\mathcal{T}}\zeta(y)\frac{\delta}{\delta\psi(y)}\\{\mathcal{A}^{-}(y|\psi)\\}dy=\frac{d}{d\epsilon}\mathcal{A}^{-}(x|\psi+\epsilon\zeta)\bigg{|}_{\epsilon=0}$

can be thought of as the functional analog of a directional derivative of
$\mathcal{A}^{-}(x|\psi)$ ‘in the direction of the field’ $\zeta$. We will now
use this linear approximation to study phenotypic clustering/evolutionary
branching using spectral methods first introduced for various specific models
by [121] and further expanded upon for some specific models of adaptive
dynamics in [122].

###  S7.2 Detecting phenotypic clusters through Fourier analysis

We begin with the linear functional Fokker-Planck equation

$\frac{\partial P}{\partial
t}(\zeta,t)=\int\limits_{\mathcal{T}}\left(-\frac{\delta}{\delta\zeta(x)}\left\\{\mathcal{D}_{\zeta}[\mathcal{A}^{-}](x)P(\zeta,t)\right\\}+\frac{1}{2}\mathcal{A}^{+}(x|\psi)\frac{\delta^{2}}{\delta\zeta(x)^{2}}\\{P(\zeta,t)\\}\right)dx$
(S81)

for describing stochastic fluctuations $\zeta$ from the deterministic solution
obtained by solving Eq. 19. Our goal is now to find a method to effectively
detect and describe evolutionary branches (modes in trait space, corresponding
to individual morphs) for this process. We will do this by measuring the
autocorrelation of the field $\phi$, a task made easier by moving to Fourier
space.

Figure S1: Schematic description of Fourier decomposition. Any density field
$\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\phi$ (shown
in red) can be decomposed into a sum of infinitely many Fourier modes
$\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{k}$
(shown in blue). In the Fourier dual space, we can look at the peaks of each
of these Fourier modes: The magnitude of
$\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{k}$
tells us how much it contributes to the actual function of interest
$\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\phi$. The
autocorrelation of
$\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\phi$ in the
trait space can be measured using the power spectral density in the Fourier
dual space.

A convenient theorem due to Weiner and Khinchin relates the autocorrelation of
a probability distribution to its power spectral density via Fourier
transformation. I will thus restrict myself to cases in which we can express
our focal function $\phi$ in terms of the Fourier basis
$\\{e^{ikx}\\}_{k\in\mathbb{Z}}$ (Figure S1). For example, if $\mathcal{T}$ is
an interval, this can be done by imposing ‘periodic boundary conditions’
(_i.e._ extending all functions from $\mathcal{T}$ to $\mathbb{R}$ in a way
that they appear periodic with period given by the length of the interval
$\mathcal{T}$). If $\mathcal{D}_{\zeta}[\mathcal{A}^{-}]$ is a translation-
invariant linear operator, then $\exp(ikx)$ acts as an eigenfunction,
significantly simplifying the calculations. I therefore assume that
$\mathcal{D}_{\zeta}[\mathcal{A}^{-}]$ takes the form:

$\mathcal{D}_{\zeta}[\mathcal{A}^{-}](x,t)=L[\zeta(x,t)]$

for a translation-invariant linear operator $L$ that only depends on $x$ and
$t$. This is not as restrictive as it initially sounds. For example, both the
Laplacian operator and the convolution operator are linear and translation
invariant. The presence of phenotypic clustering and polymorphisms can be
analyzed by examining the power spectrum of $P(\zeta,t)$ over the trait space,
which is precisely what we will do.

Let the Fourier basis representations $\zeta$, and $\mathcal{A}^{+}(x|\psi)$
be given by:

$\displaystyle\zeta(x,t)$
$\displaystyle=\sum\limits_{k=-\infty}^{\infty}e^{ikx}\zeta_{k}(t)\ \ ;\ \
\zeta_{k}(t)=\int\limits_{\mathcal{T}}\zeta(x,t)e^{-ikx}dx$ (S82)
$\displaystyle\mathcal{A}^{+}(x|\psi)$
$\displaystyle=\sum\limits_{k=-\infty}^{\infty}e^{ikx}A_{k}(t)\ \ ;\ \
A_{k}(t)=\int\limits_{\mathcal{T}}\mathcal{A}^{+}(x|\psi)e^{-ikx}dx$

In this case, the functional derivative operator obeys:

$\frac{\delta}{\delta\zeta(x)}=\sum\limits_{k=-\infty}^{\infty}e^{-ikx}\frac{\partial}{\partial\zeta_{k}}$
(S83)

and since $L$ is linear and translation-invariant, we also have the
relation222This is because $\exp(ikx)$ acts as an eigenfunction for
translation invariant linear operators, and therefore, for any function
$\varphi=\sum\varphi_{k}\exp(ikx)$, we have the relation
$L[\varphi]=L[\sum\varphi_{k}\exp(ikx)]=\sum\varphi_{k}L[\exp(ikx)]=\sum\varphi_{k}L_{k}\exp(ikx)$,
where $L_{k}$ is the eigenvalue of $L$ associated with the eigenfunction
$\exp(ikx)$. It is helpful to draw the analogy with eigenvectors of matrices
and view $L_{k}\varphi_{k}$ as the projection of $L[\varphi]$ along the $k$th
eigenvector $e_{k}=\exp(ikx)$.:

$L[\zeta]=\sum\limits_{k=-\infty}^{\infty}L_{k}\zeta_{k}e^{ikx}$ (S84)

where

$L_{k}=e^{-ikx}L[e^{ikx}]$

Lastly, by definition of Fourier modes, we have, for any differentiable real
function $F$ and any fixed time $t>0$:

$\frac{\partial}{\partial\zeta_{j}(t)}F(\zeta_{i}(t))=\delta_{i,j}F^{\prime}(\zeta_{j}(t))$
(S85)

where $\delta_{i,j}$ is the Kronecker delta, defined as

$\delta_{i,j}=\begin{cases}1&i=j\\\ 0&i\neq j\end{cases}$

Using Eq. S82, Eq. S83, and Eq. S84 in Eq. S81, we get, for the first term of
the RHS:

$\displaystyle-\int\limits_{\mathcal{T}}\frac{\delta}{\delta\zeta(x)}\left\\{L[\zeta(x,t)]P(\zeta,t)\right\\}dx$
$\displaystyle=-\int\limits_{\mathcal{T}}\sum\limits_{k}e^{-ikx}\frac{\partial}{\partial\zeta_{k}}\\{\sum\limits_{n}e^{inx}L_{n}\zeta_{n}P\\}dx$
$\displaystyle=-\int\limits_{\mathcal{T}}\sum\limits_{k}\sum\limits_{n}e^{-i(k-n)x}\frac{\partial}{\partial\zeta_{k}}\\{L_{n}\zeta_{n}P\\}dx$
$\displaystyle=-2\pi\sum\limits_{k}L_{k}\frac{\partial}{\partial\zeta_{k}}\\{\zeta_{k}P\\}$
(S86)

and for the second:

$\displaystyle\int\limits_{\mathcal{T}}\sum\limits_{k}e^{ikx}A_{k}\left(\sum\limits_{m}\sum\limits_{n}e^{-i(m+n)x}\frac{\partial}{\partial\zeta_{m}}\frac{\partial}{\partial\zeta_{n}}P\right)dx$
$\displaystyle=\int\limits_{\mathcal{T}}\sum\limits_{k}\sum\limits_{m}\sum\limits_{n}e^{i(k-m-n)x}A_{k}\frac{\partial}{\partial\zeta_{m}}\frac{\partial}{\partial\zeta_{n}}\\{P\\}dx$
$\displaystyle=2\pi\sum\limits_{m}\sum\limits_{n}A_{m+n}\frac{\partial}{\partial\zeta_{m}}\frac{\partial}{\partial\zeta_{n}}\\{P\\}$
(S87)

Substituting Eq. S86 and Eq. S87 into Eq. S81, we see that the Fokker-Planck
equation Eq. S81 in Fourier space reads:

$\frac{\partial P}{\partial
t}=-2\pi\sum\limits_{k}L_{k}\frac{\partial}{\partial\zeta_{k}}\\{\zeta_{k}P\\}+\pi\sum\limits_{m}\sum\limits_{n}A_{m+n}\frac{\partial}{\partial\zeta_{m}}\frac{\partial}{\partial\zeta_{n}}\\{P\\}$
(S88)

It is important to remember that since $\zeta(x,t)$ is a stochastic process,
$\zeta_{i}$ is really a stochastic process and thus $\zeta_{i}(t)$ is actually
shorthand for the random variable $(\zeta_{i})_{t}(\omega)$, where $\omega$ is
a sample path in the Fourier dual of our original probability space.
Multiplying both sides of Eq. S88 by $\zeta_{r}$ and integrating over the
probability space to obtain expectation values, we see that

$\displaystyle\frac{d}{dt}\mathbb{E}[\zeta_{r}]$
$\displaystyle=-2\pi\sum\limits_{k}\int\zeta_{r}L_{k}\frac{\partial}{\partial\zeta_{k}}\\{\zeta_{k}P\\}d\omega+\pi\sum\limits_{m}\sum\limits_{n}A_{m+n}\int\zeta_{r}\frac{\partial}{\partial\zeta_{m}}\frac{\partial}{\partial\zeta_{n}}(P)d\omega$
(S89)

We will evaluate the terms on the RHS of Eq. S89 using integration by parts.
Recall that for any two functions $u$ and $v$ defined on a domain $\Omega$,
the general formula for integration by parts is given by:

$\int\limits_{\Omega}\frac{\partial u}{\partial
x_{i}}vd\mathbf{x}=-\int\limits_{\Omega}u\frac{\partial v}{\partial
x_{i}}d\mathbf{x}+\int\limits_{\partial\Omega}uv\gamma_{i}dS(\mathbf{x})$
(S90)

where $\partial\Omega$ is the boundary of $\Omega$, $dS$ is the surface
element of this boundary, and $\gamma_{i}$ is the $i\textsuperscript{th}$
component of the unit outward normal to the boundary. In our case, I assume
that the probability of extreme events is negligible enough that $P$ decays
rapidly near the boundaries and we can neglect the contributions of the
boundary term (second term on the RHS of Eq. S90). Thus, using integration by
parts and neglecting the boundary terms on the RHS of equation Eq. S89, we
obtain

$\displaystyle\frac{d}{dt}\mathbb{E}[\zeta_{r}]$
$\displaystyle=2\pi\sum\limits_{k}L_{k}\int\zeta_{k}\frac{\partial\zeta_{r}}{\partial\zeta_{k}}Pd\omega+\pi\sum\limits_{m}\sum\limits_{n}A_{m+n}\int\frac{\partial^{2}\zeta_{r}}{\partial\zeta_{m}\partial\zeta_{n}}Pd\omega=2\pi
L_{r}\mathbb{E}[\zeta_{r}]$ (S91)

where we have used the relation in Eq. S85 to arrive at the final expression.
Similarly, multiplying Eq. S88 by $\zeta_{r}\zeta_{s}$, integrating over the
probability space and using integration by parts, we get:

$\displaystyle\frac{d}{dt}\mathbb{E}[\zeta_{r}\zeta_{s}]$
$\displaystyle=2\pi\sum\limits_{k}L_{k}\int\zeta_{k}P\frac{\partial}{\partial\zeta_{k}}\\{\zeta_{r}\zeta_{s}\\}d\omega+\pi\sum\limits_{m}\sum\limits_{n}A_{m+n}\int
P\frac{\partial}{\partial\zeta_{m}}\frac{\partial}{\partial\zeta_{n}}\\{\zeta_{r}\zeta_{s}\\}d\omega$
$\displaystyle=2\pi(L_{r}+L_{s})\mathbb{E}[\zeta_{r}\zeta_{s}]+\pi(A_{2r}+A_{2s})$
(S92)

At the stationary state, the LHS must be zero by definition, and we must
therefore have, for every $r,s\in\mathbb{Z}$,:

$\mathbb{E}[\zeta_{r}\zeta_{s}]=-\frac{A_{2r}+A_{2s}}{2(L_{r}+L_{s})}$ (S93)

Recall now that the Fourier modes of any real function $\varphi$ must satisfy
$\varphi_{-r}=\overline{\varphi}_{r}$. Since $\zeta$, $A$ and $L$ are all
real, we can substitute $s=-r$ in equation Eq. S93 to obtain the
autocovariance relation:

$\mathbb{E}[|\zeta_{r}|^{2}]=-\frac{\mathrm{Re}(A_{2r})}{2\mathrm{Re}(L_{r})}$
(S94)

The presence of phenotypic clustering can be detected using the ‘spatial
covariance’ of our original process $\phi$, defined as [121]:

$\Xi[x]=m(\mathcal{T})\int\limits_{\mathcal{T}}\mathbb{E}[\phi_{\infty}(x)\phi_{\infty}(y-x)]dy$
(S95)

where $\phi_{\infty}$ is the stationary state distribution of
$\\{\phi_{t}\\}_{t}$ and $m$ is the Lebesgue measure. We can use a spatial
analogue of the Wiener-Khinchin theorem to calculate [121, 122]:

$\Xi[x]=m(\mathcal{T})\bigg{[}\underbrace{\int\limits_{\mathcal{T}}\psi_{\infty}(x)\psi_{\infty}(y-x)dy}_{\text{Infinite
population
prediction}}+\underbrace{\vphantom{\int\limits_{\mathcal{T}}}\frac{1}{K}\sum\limits_{r=-\infty}^{\infty}\mathbb{E}[|\zeta_{r}|^{2}]e^{irx}}_{\text{Finite
population corrections}}\bigg{]}$ (S96)

where the expectations in the second term are for the stationary state. A flat
$\Xi[x]$ indicates that there are no clusters, and peaks indicate the presence
of clusters. Notice that we generically expect stochastic fluctuations to
induce some finite population corrections to the spatial covariance.

##  S8 An example: Asexual resource competition

In this section, I present an example of a model of resource competition to
illustrate the analytical pipeline outlined in this paper. Let us imagine a
population of asexual individuals bearing some quantitative trait $x$ taking
values in $\mathbb{R}$. I assume the particular birth and death rate
functionals

$\displaystyle b(x|\nu)$ $\displaystyle=\int M(x,y)\nu(y,t)dy\ \ ;\ \
M(x,y)=\exp\left(\frac{-(x-y)^{2}}{\sigma^{2}_{m}}\right)$ (S97)
$\displaystyle d(x|\nu)$
$\displaystyle=\frac{\nu(x,t)}{K}\int\alpha(x,y)\nu(y,t)dy\ \ ;\ \
\alpha(x,y)=\exp\left(\frac{-(x-y)^{2}}{\sigma^{2}_{\alpha}}\right)$

where the integrals are over the entire real line.

The interpretation is as follows: Individuals give birth to offspring with
mutations. The effect of the mutations is parametrized by the ‘mutation
kernel’ $M$. $M(x,y)$ represents the per-capita rate at which type $y$
individuals give birth to type $x$ individuals. I assume $M$ has the
functional form of a Normal distribution centered around the focal trait value
and with a variance of $\sigma^{2}_{m}$. Thus, birth is most often with no
mutation, and more extreme mutations are less likely. The death rate
functional incorporates the effect of death due to resource competition. The
effects of competition are parametrized by a competition kernel $\alpha$. The
quantity $\alpha(x,y)$ represents the per-capita additional death rate of type
$x$ individuals due to competition with type $y$ individuals. I have assumed
that the strength of competition effects takes the form of a Normal
distribution about the focal trait value, with variance
$\sigma^{2}_{\alpha}>0$. Thus, individuals experience less competition from
individuals who are at a greater phenotypic distance away from them in trait
space. The death rate term also incorporates a carrying capacity $K>0$. This
functional form for the death rate can be seen as the continuous analog of the
death rate term that arises in Lotka-Volterra competition
($x_{i}\sum_{j}\alpha_{ij}x_{j}/K$).

Switching from population numbers $\nu$ to the population density field
$\phi=\nu/K$ via the variable transformation presented in Eq. S2, we obtain
the new stochastic field equations paramterized by the scaled birth and death
rate functionals:

$\displaystyle b(x|\phi)$ $\displaystyle=\int M(x,y)\phi(y,t)dy$ (S98)
$\displaystyle d(x|\phi)$ $\displaystyle=\phi(x,t)\int\alpha(x,y)\phi(y,t)dy$

Assuming the mutational variance $\sigma^{2}_{m}$ is small, we can approximate
the birth rate functional as

$\int
e^{-(x-y)^{2}/{\sigma^{2}_{m}}}\phi(y,t)dy=\phi(x,t)+\frac{\sigma^{2}_{m}}{2}\frac{\partial^{2}\phi}{\partial
x^{2}}+\ldots$ (S99)

Discarding higher order terms, we can thus write the birth and death rates
presented in Eq. S98 as

$\displaystyle b(x|\phi)$
$\displaystyle=\phi(x,t)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}\phi$ (S100)
$\displaystyle d(x|\phi)$ $\displaystyle=\phi(x,t)\int\alpha(x,y)\phi(y,t)dy$

where we have used the notation $\nabla^{2}_{x}=\partial^{2}/\partial x^{2}$
for convenience. We can now compare terms with the definitions of fitness,
turnover, and mutational effects to identify:

$\displaystyle w(x|\phi)$ $\displaystyle=1-\int\alpha(x,y)\phi(y,t)dy$ (S101)
$\displaystyle\tau(x|\phi)$ $\displaystyle=1+\int\alpha(x,y)\phi(y,t)dy$
(S102) $\displaystyle\mu$ $\displaystyle=\sigma^{2}_{m}/2$ (S103)
$\displaystyle Q(x|\phi)$ $\displaystyle=\nabla^{2}_{x}\phi$ (S104)

These quantities are sufficient to derive all the general equations presented
in this paper. For conciseness, I only present the resultant equations for the
density field (Eq. 17) and the stochastic gradient equation (Eq. 30) below.

###  S8.1 The population density field

From Eq. 17, we see that the population density field obeys the SPDE

$\begin{split}\frac{\partial\phi}{\partial
t}(x,t)&=\phi(x,t)\left(1-\int\alpha(x,y)\phi(y,t)dy\right)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}\phi\\\
\hphantom{=}&+\frac{1}{\sqrt{K}}\sqrt{\phi(x,t)\left(1+\int\alpha(x,y)\phi(y,t)dy\right)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}\phi\hphantom{+}}\dot{W}(x,t)\end{split}$
(S105)

Equation S105 is a stochastic version of the quantitative logistic equation
that includes the effects of noise-induced selection and drift. This can be
seen by examining the infinite population limit.

#### The infinite population limit

In the infinite population limit, we see from Eq. 19 that the density field
obeys the PDE

$\frac{\partial\phi}{\partial
t}(x,t)=\phi(x,t)\left(1-\int\alpha(x,y)\phi(y,t)dy\right)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}\phi$
(S106)

Equation S106 is the so-called ‘quantitative logistic equation’ and has been
used to model asexual resource competition [106].

###  S8.2 The stochastic gradient equation

Let us first calculate the relevant fitness and turnover functions. From Eq.
S101, we obtain:

$w(x|N_{K}\delta_{y})=1-\int\alpha(x,z)N_{K}\delta_{y}dz=1-N_{K}\alpha(x,y)$
(S107)

and similarly, from Eq. S102, we obtain

$\tau(x|N_{K}\delta_{y})=1+\int\alpha(x,z)N_{K}\delta_{y}dz=1+N_{K}\alpha(x,y)$
(S108)

Substituting Eq. S107 and Eq. S108 into Eq. S74, we thus obtain the stochastic
gradient equation:

$\frac{dy}{dt}=\sigma^{2}_{x}(t)\frac{\partial G(x;y)}{\partial
x}\bigg{|}_{x=y}+\frac{dW_{y}}{dt}$ (S109)

where

$G(x;y)=1-\frac{1}{KN_{K}(t)}-\left(\frac{1+KN_{K}(t)}{K}\right)\alpha(x,y)$
(S110)

and

$\frac{dW_{y}}{dt}=\int(x-y(t))\sqrt{1+N_{K}(t)\alpha(x,y)}\dot{W}(x,t)dx$
(S111)

Equation S109 is a gradient equation that incorporates the effects of noise-
induced selection and genetic drift, and thus can be used to study evolution
in finite populations. The infinite population limit recovers a well-known
model from adaptive dynamics, as I show below.

#### The infinite population limit

In the infinite population limit, we see from Eq. S74 and Eq. S107 that we
obtain the gradient equation

$\frac{dy}{dt}=-\sigma^{2}_{x}(t)\frac{\partial\alpha(x,y)}{\partial
x}\bigg{|}_{x=y}$ (S112)

Equation S112 is the asexual model of sympatric speciation due to resource
competition presented in Chapter 3 of [106].

##  S9  A stochastic Fisher-KPP equation from an individual-based model

Consider the birth and death rate functionals given by S97 with the
competition kernel $\alpha(x,y)=\delta_{x}$. Thus, we can write

$\displaystyle b(x|\nu)$ $\displaystyle=\int M(x,y)\nu(y,t)dy\ \ ;\ \
M(x,y)=\exp\left(\frac{-(x-y)^{2}}{\sigma^{2}_{m}}\right)$ (S113)
$\displaystyle d(x|\nu)$
$\displaystyle=\frac{\nu(x,t)}{K}\int\delta_{x}\nu(y,t)dy=\frac{1}{K}\nu^{2}(x,t)$

Biologically, $\alpha(x,y)=\delta_{x}$ means that individuals only compete
with other individuals that have the exact same phenotype as them, and can be
thought of as taking the limit of infinitely narrow niche width in the
competition model introduced in the previous section. For instance, such
dynamics could be a reasonable model if the trait in question is location in
physical space. Using the same approximation for the birth rate as in the
previous section, we arrive at the birth and death rates in density space

$\displaystyle b(x|\phi)$
$\displaystyle=\phi(x,t)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}\phi$ (S114)
$\displaystyle d(x|\phi)$ $\displaystyle=\phi^{2}(x,t)$

and thus

$\displaystyle w(x|\phi)$ $\displaystyle=1-\phi(x,t)$ (S115)
$\displaystyle\tau(x|\phi)$ $\displaystyle=1+\phi(x,t)$ (S116)
$\displaystyle\mu$ $\displaystyle=\sigma^{2}_{m}/2$ (S117) $\displaystyle
Q(x|\phi)$ $\displaystyle=\nabla^{2}_{x}\phi$ (S118)

We can now subtitute these functions into the population level equations for
population densities (abundances).

###  S9.1 The population density field

From Eq. 17, we see that the population density field obeys the SPDE

$\begin{split}\frac{\partial\phi}{\partial
t}(x,t)&=\phi\left(1-\phi\right)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}\phi+\frac{1}{\sqrt{K}}\sqrt{\phi\left(1+\phi\right)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}\phi\hphantom{+}}\dot{W}(x,t)\end{split}$
(S119)

Equation S119 is similar to a stochastic analog of the Fisher-KPP equation on
$\mathbb{R}$, as will become clear upon taking the infinite population limit
below. However, note that Eq. S119 is _not_ the SPDE that is usually referred
to as the ‘stochastic Fisher-KPP’ equation [107, 96]. The usual ‘stochastic
Fisher-KPP’ equation is instead the SPDE [107, 96]

$\begin{split}\frac{\partial u}{\partial
t}(x,t)&=u\left(1-u\right)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}u+\epsilon\sqrt{u(1-u)}\dot{W}(x,t)\end{split}$
(S120)

The discrepancy arises because the SPDE (S120) in fact does _not_ correspond
to the mesoscopic limit of a (spatially-extended/infinite dimensional) birth-
death process [107] but is instead written down directly in analogy to the
discrete trait Wright-Fisher diffusion

$dp=p(1-p)dt+\sqrt{p(1-p)}dW_{t}$ (S121)

While Eq. S120 can be connected with birth-death processes via duality [107]
and can be derived as a scaling limit of either a contact process [115] or a
spatial $\Lambda$-Fleming-Viot process [96], I show here that it does not
arise as the expected mesoscopic limit of a measure-valued _birth-death_
process. Instead, the mesoscopic limit I obtain is Eq.S119, which is also the
equation expected from the more rigorous measure-theoretic approach to
studying the processes I consider in this paper (see section 4.2 in [99] with
$\eta=1$).

#### The infinite population limit

In the infinite population limit, we see from Eq. 19 that the density field
obeys the PDE

$\frac{\partial\phi}{\partial
t}=\phi\left(1-\phi\right)+\frac{\sigma^{2}_{m}}{2}\nabla^{2}_{x}\phi$ (S122)

This is the famous Fisher-KPP equation and is widely used to model spatial
evolution [119, 97, 96].


# EvoCut : A new Generalization of Albert-Barabási Model for Evolution of
Complex Networks

Shailesh Kumar Jaiswal

Nabajyoti Medhi

Manjish Pal

National Institute of Technology  Meghalaya. 793003<EMAIL_ADDRESS>Tezpur  University  Assam - 784028<EMAIL_ADDRESS>National Institute of
Technology  Meghalaya. 793003<EMAIL_ADDRESS>Mridul Sahu

Prashant Sahu

Amal Dev Sarma

National Institute of Technology  Meghalaya. 793003<EMAIL_ADDRESS>National Institute of Technology  Meghalaya. 793003<EMAIL_ADDRESS>National Institute of Technology  Meghalaya. 793003<EMAIL_ADDRESS>
(24 April 2018)

###### Abstract

With the evolution of social networks, the network structure shows dynamic
nature in which nodes and edges appear as well as disappear for various
reasons. The role of a node in the network is presented as the number of
interactions it has with the other nodes. For this purpose a network is
modeled as a graph where nodes represent network members and edges represent a
relationship among them. Several models for evolution of social networks has
been proposed till date, most widely accepted being the Barabási-Albert [1]
model that is based on _preferential attachment_ of nodes according to the
degree distribution. This model leads to generation of graphs that are called
_Scale Free_ and the degree distribution of such graphs follow the _power
law_. Several generalizations of this model has also been proposed. In this
paper we present a new generalization of the model and attempt to bring out
its implications in real life.

###### keywords:

Social Media Networks, Evolution of Social Networks, Scale Free Graphs,
Barabási-Albert model.

††conference:  Lyon, France††terms: WWW, Social media networks, Online social
networks.

##  1 Introduction

Our lives are surrounded with several complex phenomena which are at times
very difficult to explain or understand. These phenomena usually encompass a
lot of spheres of our lives. One such phenomena is how certain biological,
geological, physical, astronomical, financial and social systems show a very
peculiar similarity that is they all exhibit a so called _Scale Free Property_
[4]. This property says that certain features of these systems follow a common
pattern. For example in case of geology, most of the earthquakes that occur on
this planet are nominal whereas few are gigantic, in case of financial systems
most people on earth have small incomes and a few have a monumental incomes,
in case of social systems, there are few celebrities who are extremely popular
whereas most of the other are not and this list of examples continues on and
on. As it turns out this kind of a Scale free pattern has been known to
scientists and researchers for long. Herbert Simon [3] showed the existence of
such properties way back in the 1950s. Although this property is very common,
present literature exhibits that there is not much understanding of as to why
and how it occurs. In recent years with the advent of new technologies like
the internet, WWW, social media etc., a new stream of research has come up
which seeks to explain this phenomenon in several complex systems by employing
mathematical models. Here the key is to find the existence of a graph
theoretic structure that is present in these systems and show that the Scale
Free Property is present in that graph. Researchers have shown that the
presence of this property means that the corresponding graph follows the
_power law_ degree distribution. The most well known among such models is the
model of Barab$\acute{a}$si-Albert that describes a model based on
_preferential attachment_. This model ensures that as this graph evolves by
the addition of new nodes, each new node gets attached to a set of $m$
existing nodes where the attaching probability of the new node is proportional
to the degree of that node already present in the graph. This model although
considered to be fairly satisfactory, is not very realistic at times and has
several drawbacks. Due to this several generalizations of this model has been
proposed that try to improvise on the properties of the graph generated. In
this paper we propose a new model that is a generalization of the AB model and
is based on the _cuts_ in the graph. Our model is very novel and has no
apparent links with the already existing generalization of the Barabási-Albert
(BA) model. In this paper, we introduce our model, which we have named as the
_EvoCut_ model, and describe its properties. We further bring out how this
model is more realistic than the already existing models.

###  1.1 The Scale Free Property

As mentioned in the previous section, the scale free property [4] basically
means that the occurrence of very high is small and there is an abundance of
small. In the context of networks and graphs this means that there are very
few nodes with high degree and there are a lot of nodes with small degree. A
more technical way of putting it is to say that the degree distribution
follows the _power law_ i.e. $p(k)\propto k^{-\gamma}$ where $\gamma\in(2,3]$;
$p(k)$ is the number of nodes with degree $k$ divided by $2m$ (which is the
total degree of nodes). This alternatively means that on log scale the plot of
$p(k)$ vs. $k$ is linear.

##  2 Prior Work on Generalization of Albert-Barabási model

In this section we present some prior models that are based on the
generalization of the Barab$\acute{a}$si-Albert model.

###  2.1 Dorogovtsev-Mendes-Samukhin (DMS)
Model

Dorogovtsev, Mendes and Samukhin present a model [5] that generalizes the BA
model in the following sense, they incorporate the initial state of the nodes
in the network and the degree ($k$) is dependent on the initial state ($s$)
and time ($t$). They come up with a dynamics that achieves power law when the
degree is very large i.e. $p(k)\propto k^{-\gamma}$ when $k\rightarrow\infty$
and $\gamma\in(2,\infty)$. The average connectivity
$k(s,t)\propto(s/t)^{-\beta}$ as $(s/t)\rightarrow 0$. This generalization
makes the model very unrealistic because $\gamma\in(2,\infty)$ whereas in
reality $\gamma\in(2,3]$, as supported by AB model.

###  2.2 Antal-Krapivsky-Redner (AKR) Model

The AKR model [6] generalizes the BA model by considering a graph in which the
links are associated with a notion of friendship (+1) or enmity (-1) and the
notion of a node is replaces by a collection of three nodes linked with either
(+1) or (-1) links, which is called a triad. The model defines a notion of a
balance which is the product of all the values of the links in a particular
triad. If this product is equal to 1 then it is called _balanced_ otherwise it
is called _imbalanced_. In this model, the authors define a dynamics that
tries to maintain the balance of the triads and describes how the degree of a
link (which is the number of triads in which it is involved) changes with time
as the network is allowed to evolve according to the dynamics. They derive a
set of differential equations that governs this rate and solve them to get a
distribution same as the power law.

###  2.3 Sole-Pastor-Satorras-Smith-Kepler (SPSK) Model

The Sol$\acute{e}$, Pastor-Satorras, Smith, Kepler (SPSK) model [7] uses 3
mechanisms _duplicate_ , _divergence_ and _mutate_ for the process of
evolution. Using the operation _Duplicate_ one can copy a randomly selected
node along with its connections, using _Divergence_ operation one can delete
some connections made after the duplicate operation and the _Mutate_ operation
allows us to add connections once the duplicate operation is applied. It has
been shown by them that this generalization produces the _power law_ but leads
to unrealistic assortativity and clustering. This is a generalization of the
BA model because it allows more possible ways of connection of a new node with
the existing graph.

It is customary to note that there have been some generalizations of
preferential attachment model that are based on the addition and deletion of
nodes [8, 9].

##  3 Our Contribution

In this paper we present a new model that generalizes the BA model. As
compared to the previous models our model is a natural generalization of the
BA model and is more intuitive than the previous generalizations. It uses the
combinatorial properties of the given graph and does not create unrealistic
links as in the case of SPSK model and DMS model. Our model has two variants,
one of which generates scale free networks (as in BA model) and the other one
gives rise to a family of graphs in which the degree distribution follows the
stretched exponential distribution rather than the power law, which is similar
to the SPSK model. In the following sections we describe in detail our
proposed model and also discuss the properties of the two variants of the
model.

##  4 Preliminaries

We model the network as an undirected graph in which nodes represent the
network members and an undirected edge represents a relationship between them.
Initially the graph is considered to have $n_{0}$ number of nodes with a few
links $m_{0}$. Our evolution process crucially uses the notion of cuts which
is an important combinatorial property of a graph. Given a graph $G=(V,E)$ a
subset of vertices $S$ and $\bar{S}=V\backslash S$, a _Cut_ is defined as the
set $E(S,\bar{S})$ $=\\{(a,b)|a\in S\mbox{ and }b\in\bar{S}\\}$. Our evolution
process uses the size of the cut crucially making it substantially different
from all the prior models that do not use the notion of cut in any way. We
also define a $k$-neighborhood of a particular vertex $v$ as the number of
nodes of $G$ which are with in a distance $k$ from the node $v$. Thus
$B(v,k)=\\{u|d(v,u)\leq k\\}$. We also define the set $B^{\prime}(v,k)$ as the
set

$\\{u\in B(v,k)|(u,v)\in E(B(v,k),\overline{B(v,k)});v\in\overline{B(v,k)}\\}$

##  5 The EvoCut Model

In this section we introduce our model and describe in detail its two
variants. In the first variant a new node is attached to a node which has the
maximum pulling power which is based on the size of the $k$-neighborhood of
that node whereas in the second variant the node is attached to a randomly
chosen node which is on the boundary of the $k$-neighborhood set of the node
with maximum pulling power.

$G_{0}=(V_{0},E_{0})$ ; $t=0$; $|V_{0}|=n_{0}$ and $E_{0}=m_{0}$, $k=k_{o}$,
$Y=0$, $m=0$;

while  _( $n_{t}\leq N$)_ do

Let $v_{t}$ be the new node at time $t$;

for  _$v\in V_{t}$_ do

Compute $x=|E(B(v,k),\overline{B(v,k)}|$;

$Y=Y+x$;

end for

Compute $x_{v}=\mbox{max}_{v}\\{\frac{|E(B(v,k),\overline{B(v,k)}|}{Y}$};

Compute $v^{\prime}=\mbox{argmax}_{v}\\{x_{v}\\}$;

$E_{t}=E_{t}\bigcup(v_{t},v^{\prime});V_{t}=V_{t}\bigcup v_{t}$;

$t=t+1$;

$Y=0$;

end while

Algorithm 1 Model A

$G_{0}=(V_{0},E_{0})$ ; $t=0$; $|V_{0}|=n_{0}$ and $E_{0}=m_{0}$, $k=k_{o}$,
$Y=0$, $m=0$;

while  _( $n_{t}\leq N$)_ do

Let $v_{t}$ be the new node at time $t$;

for  _$v\in V_{t}$_ do

Compute $x=|E(B(v,k),\overline{B(v,k)}|$;

$Y=Y+x$;

end for

Compute $x_{v}=\mbox{max}_{v}\\{\frac{|E(B(v,k),\overline{B(v,k)}|}{Y}$};

Select a node $v^{\prime}$ randomly uniformly from the set $B^{\prime}(v,k)$;

$E_{t}=E_{t}\bigcup(v_{t},v^{\prime});V_{t}=V_{t}\bigcup v_{t}$;

$t=t+1$;

$Y=0$;

end while

Algorithm 2 Model B

##  6 Properties of the EvoCut Model

In this section we describe in detail the properties and the reasoning behind
coming up with these models. In this model we allow the nodes to be attached
one by one and the pulling power of a particular node is defined by the size
of the cut $E(B(v,k),\overline{B(v,k)}$ which the number of edges in the
$k$-neighborhood of the node $v$ where $k$ is a parameter between $[0-(n-1)]$.
The pulling power is same for both the models $A$ and $B$. The other features
of the models is being mentioned in the following subsections.

###  6.1 Model A - Deterministic Case

In this model a particular node is attached to that node which has the maximum
pulling power based on a given value of $k$. This model is thus fully
_deterministic_ i.e. doesn’t use any source of randomness. The process starts
with an initial graph $G_{0}=(V_{0},E_{0})$ and goes until the number of nodes
in the graph is less than $N$. Every time when a new node arrives, the for
loop computes the normalizing factor

$\displaystyle\sum_{v}|E(B(v,k),\overline{B(v,k)}|$ (1)

and the later part of the model computes the node $v$ that maximizes the ratio

$\displaystyle\frac{B(v,k)}{\sum_{u}|E(B(u,k),\overline{B(u,k)}|}$ (2)

and attaches the new node with this node.

###  6.2 Model B - Randomized Case

In this model a particular node is attached to a randomly chosen node on the
boundary of the set $B(v,k)$. The boundary is defined as the nodes in $B(v,k)$
which are incident to the edges in the set $E(B(v,k),\overline{B(v,k)}$. This
model thus uses randomness in choosing the node to which the incoming node
should be added. As in the previous case the model first computes the node $v$
that maximizes the factor

$\displaystyle\frac{B(v,k)}{\sum_{u}|E(B(u,k),\overline{B(u,k)}|}$ (3)

and then randomly chooses a node from the set $B^{\prime}(v,k)$ and attach the
new node with that node.

###  6.3 As a Generalization of BA Model

One can observe that in both the models mentioned above if we fix the
parameter $k$ as 0, then we get the (BA) model. This is simply because when
$k=0$, $B(v,k)=v$ and the set $E(B(v,k),\overline{B(v,k)}$ is equal to the
degree of $v$ and in that case the attaching probability distribution is
exactly equal to the degree distribution of the graph at time $t$ which is
same as the BA model. Thus, in our model, only the special case of $k=0$
results in the BA model.

##  7 Experimental Results and Analysis

In this section we present the degree distributions that we obtain as we let
our model to evolve on a set of nodes with certain initial condition.

###  7.1 Analysis of Model A

The following are the plots that we obtain once the degree distributions are
generated for _Model A_.

According to the plots that are generated we can infer the following:

  * •

For small values of $k$ this model gives rise to a scale free distribution
when $k$ is even and for larger values of $k$ the degree distribution follows
the _stretched exponential distribution_ in the log-scale. This in turn
implies that in the normal scale the plot ensures that the number of nodes
with high degree is fairly large as compared to the _scale free distribution_
which implies that for large values of $k$ there are several nodes with high
degree.

  * •

The intuitive reasoning for the aforementioned result is that for small values
of $k$ this model behaves basically similar to the BA model. Whereas for
larger values of $k$ we can observe that once a new node gets attached to the
already existing node the pulling power of the already existing node does not
increase and the pulling power of the $k^{th}$-neighbor increases. This
phenomenon creates some sort of an oscillation on the increase in the degree
of the certain number of nodes which implies that this set of nodes experience
enhancement of degree simultaneously. This explains the _stretched exponential
distribution_ in the log-scale for larger values of $k$.

###  7.2 Analysis of Model B

The following are the plots that we obtain once the degree distributions of
_Model B_ are generated and based on these we infer the following

  * •

As said earlier in this model the incoming node is attached to that node which
is on the boundary of the $k$-neighborhood of the node that maximizes the
objective function.

  * •

From the plots, it is clear that for both small and large values of $k$ we get
a _power law_ distribution. This observation can be explained intuitively as
follows: since the incoming node $x$ is getting attached to a node which is on
the boundary of the $k$-neighborhood of a node $y$, the pulling power of node
$y$ increases with time because it is defined as the number of edges on the
boundary of the $k$-neighborhood of $y$ and hence in further iterations the
node $y$ gets enhancement in its power.

  * •

This makes the model a generalization of the BA model which is quite similar
to it.

##  8 Comparison with Prior Models

In this section we compare our models with the prior generalizations of the BA
model and also discuss how close to reality these models can be considered.

###  8.1 Comparison with DMS Model

As mentioned before, the DMS model starts with an initial state and the degree
of the node is dependent on the initial state of the evolution process, a
feature which our models also possess. The DMS model achieves the power law
for large values of the degree whereas in our case the Model B achieves power
law even for small values of the degree.

###  8.2 Comparison with AKR Model

In the AKR model the authors come up with a notion of a balance of the values
present in a triad which is a collection of three nodes and the directed links
among them. The dynamics of the model tries to maintain the balance of the
triads. This model, although a generalization of the BA model uses only
_local_ modifications whereas our models are _global_ in the sense that each
node looks around a $k$-neighborhood where $k$ can be fairly large. This
allows our model to look at the global influence of the nodes which is not
present in the AKR model.

###  8.3 Comparison with SPSK Model

In the SPSK Model, the evolution of the graph is based on some operations done
on the original graph called _duplicate_ , _divergence_ and _mutate_. This
make the model very restrictive because a node can only get attached to the
existing graph by performing a duplicate operation which is basically the
replication of the connection of an already existing node. Our models are not
restrictive in that sense and a new node is not forced to follow the topology
of an already existing node and is attached purely based on the pulling power
of a particular node. Despite the difference in nature of SPSK and our models
we observe that the degree distribution of this model and Model A proposed by
us is quite similar and follows the stretched exponential distribution.

##  9 Real Life Implications of EvoCut

One of the important questions regarding the models of evolution of complex
networks is how realistic are they. The BA model although satisfactory is not
considered to be very realistic. Thus there is definitely a need to come up
with models which not only give rise to the power law distribution but are
also realistic. In this section we justify how our models can explain the
realities of certain complex networks. In the case of Model A, as we had
observed that this model leads to a degree distribution [10] in which there
are lots of nodes with large degrees. Although this doesn’t follow the scale
free property but explains the nature of some real political networks. It is
well known that [11] the stretched exponential distribution describe very well
the distributions of radio and light emissions from galaxies, of country
population sizes, of daily Forex US-Mark and Franc-Mark price variations, of
Vostok temperature variations and of citations of the most cited physicists in
the world. This type of distribution may not be considered to be explainable
by the BA model and its generalizations which tend to result in a power law
distribution, but can be explained by Model A.
In the case of Model B we notice that since the model looks at a
$k$-neighborhood of a particular node to compute its pulling power this allows
us to look at a larger influence of a node as compared to the standard BA
model. In fact in realistic scenarios one can conceive of several situations
in which when a new node comes it doesn’t get attached to the node with
maximum pulling power but gets attached to a node which is at a certain
distance from the node. Consider the example of celebrities in real life
networks, when a person considers himself as a follower of a particular
celebrity then he gets himself linked to another follower of the same
celebrity rather than directly getting linked with the celebrity, and this
follower might as well not be one of the closest to the celebrity i.e. one
gets attached to a _fan club_ rather than the celebrity himself.
We can also consider the example of any hierarchical ecosystem [12] present in
an organization like judiciary, police, banking system or a corporation in
which any customer is not allowed to connect with the highest authority in
that system but gets connected to a node that is lowest in that hierarchy.
This notion is being captured in the case of Model B, in the sense that the
$k$-neighborhood of a node basically creates a hierarchy around that node
considering that node as most powerful. As the distance from the node
increases we get down in the hierarchy, and hence when a new node comes under
the influence of the $k$-neighborhood of a particular node at the $k^{th}$
level of the hierarchy.

Consider the aforementioned figure in which the CEO represents the node whose
pulling power is being computed and the 1-neighborhood of it consists of VP of
Sales and VP of Service, the 2-neighborhood consists of Sales Managers and
Service Managers and the 3-neighborhood consists of Sales and Support. Thus
when an new customer comes it is likely to get in touch with the
3-neighborhood rather than the CEO. Thus in several scenarios Model B gives a
better picture of the evolution of the network and hence is much more
realistic than the BA model.

##  10 Conclusion and Future Work

In this paper, we have presented a new generalization of the Barabási-Albert
(BA) model that is based on the neighborhood properties of nodes in the
evolving graph. This generalization is substantially different from the
previous generalizations of the BA model and one of its variant gives rise to
the power law distribution and can model the growth of several real life
networks which don’t seem to be explainable by the BA model. The other variant
of our Model, despite the fact doesn’t give the scale free property models the
growth of some other real life networks and gives rise to a stretched
exponential distribution. We believe that our model is very powerful and
further investigation into it can make it a stronger candidate in
understanding the growth of several complex networks.
As part of future work it would interesting to get a mathematical proof of the
distribution generated by our models and also discover new interesting
variants of our model.

## References

  * [1] A. L. Barabási, Network science. _Cambridge university press_ , 2016.
  * [2] A. L. Barabási, Scale-free networks: a decade and beyond. science, 325 (5939), pp.412-413, 2009.
  * [3] H. A. Simon, A behavioral model of rational choice. The quarterly journal of economics, 69(1), 99-118, 1955.
  * [4] A. L. Barabási & E. Bonabeau, Scale-free networks. Scientific american, 288(5), 60-69, 2003.
  * [5] S. N. Dorogovtsev, J. F. F. Mendes & A. N. Samukhin, Structure of Growing Networks: Exact Solution of the Barabási–Albert’s Model, 2000.
  * [6] T.Antal, P. L. Krapivsky & S. Redner, Dynamics of social balance on networks. Physical Review E, 72(3), 036121, 2005.
  * [7] R. V. Solé, R. Pastor-Satorras, E. Smith & T. B. Kepler, A model of large-scale proteome evolution. Advances in Complex Systems, 5(01), 43-54, 2002.
  * [8] H. Bauke, C. Moore, J. B. Rouquier & D. Sherrington, Topological phase transition in a network model with preferential attachment and node removal. The European Physical Journal B, 83(4), 519-524, 2011.
  * [9] C. Moore, G. Ghoshal & M. E. Newman, Exact solutions for models of evolving networks with addition and deletion of nodes. Physical Review E, 74(3), 036121, 2006.
  * [10] L. A. Adamic, R. M. Lukose, A. R. Puniyani & B. A. Huberman, Search in power-law networks. Physical review E, 64(4), 046135, 2001.
  * [11] J. Laherrere & D. Sornette, Stretched exponential distributions in nature and economy:“fat tails” with characteristic scales. The European Physical Journal B-Condensed Matter and Complex Systems, 2(4), 525-539, 1998.
  * [12] E. Ravasz, A. L. Somera, D. A. Mongru, Z. N. Oltvai & A. L. Barabási, Hierarchical organization of modularity in metabolic networks. science, 297(5586), 1551-1555, 2002.