# COMPUTATIONAL PROBABILITY AND MATHEMATICAL MODELING

EDITED BY : José Roberto Cantú-González, F-Javier Almaguer, Javier Morales-Castillo and Pavel Solin PUBLISHED IN : Frontiers in Applied Mathematics and Statistics and Frontiers in Genetics

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-244-2 DOI 10.3389/978-2-88963-244-2

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# COMPUTATIONAL PROBABILITY AND MATHEMATICAL MODELING

Topic Editors:

José Roberto Cantú-González, Universidad Autónoma de Coahuila, Mexico F-Javier Almaguer, Universidad Autónoma de Nuevo León, Mexico Javier Morales-Castillo, Universidad Autónoma de Nuevo León, Mexico Pavel Solin, University of Nevada, Reno, United States of America

Image: Computational Probability and Mathematical Modeling. By: José Roberto Cantú-González and F-Javier Almaguer.

In the present time, two of the most important approaches to tackle complex systems are probability and stochastic processes theory. Still from an analytic perspective, modeling and solving a problem using a stochastic approach is not a trivial issue, hence, a combination of the logic of probabilistic reasoning with computational science is needed to obtain qualitatively good solutions in a reasonable time.

This eBook presents an interesting view of applications associated to fields of probability, statistics, and mathematic modeling, all of them supported by a computational context though the approach of stochasticity and simulation used in most of them. This collection contains three chapters, which bring applications in fields of biology, finance and physics, each chapter contains work(s) with specific applications. An editorial is also contained with a summarized version of each work, and each of them are widely explained in a specific section, which include a state of art to support the nature of the individual research, a methodology to solve the defined problem and the results and conclusions.

We hope the present eBook can represent a potential source of knowledge for the academic community of implicated disciplines, and an inspirational starting point of starting for scientists in the amazing world of applied mathematics and the search to solve complex problems

Citation: Cantú-González, J. R., Almaguer, F-J., Morales-Castillo, J., Solin, P., eds. (2019). Computational Probability and Mathematical Modeling. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-244-2

# Table of Contents

*04 Editorial: Computational Probability and Mathematical Modeling - A Stochastic Approach in Applied Sciences* José Roberto Cantú-González and F-Javier Almaguer

#### CHAPTER 1

#### COMPUTATIONAL PROBABILITY AND MATHEMATICAL MODELING IN BIOLOGY

*07 Output-Feedback Control of Virus Spreading in Complex Networks With Quarantine*

Luis A. Alarcón-Ramos, Roberto Bernal Jaquez and Alexander Schaum


#### CHAPTER 2

#### COMPUTATIONAL PROBABILITY AND MATHEMATICAL MODELING IN FINANCE

*35 Probabilistic Assessment of Investment Options in Honey Value Chains in Lamu County, Kenya*

Joshua Wafula, Yusuf Karimjee, Yvonne Tamba, Geoffrey Malava, Caroline Muchiri, Grace Koech, Jan De Leeuw, Josephat Nyongesa, Keith Shepherd and Eike Luedeling

*46 The Amnesiac Lookback Option: Selectively Monitored Lookback Options and Cryptocurrencies*

Ho-Chun Herbert Chang and Kevin Li

#### CHAPTER 3

#### COMPUTATIONAL PROBABILITY AND MATHEMATICAL MODELING IN PHYSICS

*61 On the Aerodynamic Forces on a Baseball, With Applications*

Gerardo J. Escalera Santos, Mario A. Aguirre-López, Orlando Díaz-Hernández, Filiberto Hueyotl-Zahuantitla, Javier Morales-Castillo and F-Javier Almaguer

# Editorial: Computational Probability and Mathematical Modeling - A Stochastic Approach in Applied Sciences

#### José Roberto Cantú-González <sup>1</sup> \* and F-Javier Almaguer <sup>2</sup> \*

<sup>1</sup> Escuela de Sistemas PMRV, Campus Acuña, Universidad Autónoma de Coahuila, Saltillo, Mexico, <sup>2</sup> Facultad de Ciencias Físico Matemáticas, Universidad Autónoma de Nuevo León, San Nicolás de los Garza, Mexico

Keywords: mathematical modeling, probability, stochastic processes, simulation, statistics

#### **Editorial on the Research Topic**

INTRODUCTION

**Computational Probability and Mathematical Modeling - A Stochastic Approach in Applied Sciences**

#### Edited and reviewed by:

Daniel Potts, Technische Universität Chemnitz, Germany

#### \*Correspondence:

José Roberto Cantú-González roberto.cantu@uadec.edu.mx F-Javier Almaguer francisco.almaguermrt@uanl.edu.mx

#### Specialty section:

This article was submitted to Mathematics of Computation and Data Science, a section of the journal Frontiers in Applied Mathematics and Statistics

> Received: 30 July 2019 Accepted: 13 August 2019 Published: 29 August 2019

#### Citation:

Cantú-González JR and Almaguer F-J (2019) Editorial: Computational Probability and Mathematical Modeling - A Stochastic Approach in Applied Sciences. Front. Appl. Math. Stat. 5:45. doi: 10.3389/fams.2019.00045 Probability and the stochastic processes theory, joint to the mathematical modeling and the respective computational support are clearly important work tools to tackle complex systems, which has been the fundament to develop the Research Topic of "Computational Probability and Mathematical Modeling – A Stochastic Approach in Applied Sciences".

In a brief analysis of the elements contained inside of the Research Topic, we find three important scopes: computational probability, mathematical modeling, and stochastics processes:

Computational probability, it is defined as the development of data structures and algorithms to automate the derivation of existing and new results in probability and statistics [1], Mathematical modeling is the translation of a specific problem from the natural sciences (experimental physics, chemistry, biology, geosciences) or the social sciences, or from technology, into a well-defined mathematical problem [2] and a stochastic process is a probability model that describes the evolution of a system evolving randomly in time [3]. In addition, advances in computer science have been very determinant for growing of science and technology, and all fields of mathematics included stochastic processes have not been bypassed by the digital revolution, because numerical calculation and computer simulation play a decisive role in present-day [4].

Probability and mathematical modeling are usually very adaptable for different proposes, and in definitive, they are the starting point to solve relevant problems [5].

Specifically probability and the stochastic processes theory can be applied in any field where random plays a preponderant role and /or when an analytic treatment is too complicated to be solved. In the other context, mathematical modeling as a descendent from applied mathematics rejects the abstractionism to support the science for understanding of the nature laws; which explains that mathematical models can attend disciplines as physics, biology, finances, economics, any engineering and a variety of fields to solve real problems though its included skills. In this perspective, both probability / stochastic processes theory as the mathematical modeling converge on the same objective: attention of real problems in complex systems.

## RESEARCH TOPIC AND THE ARTICLES CONTAINED

This Research Topic of "Computational Probability and Mathematical Modeling – A Stochastic Approach in Applied Sciences" presents a little but substantial view of studies related to interesting scenarios from the applied science in fields of genetics, biology, physics, finances, and investment options. In addition, as a core and distinguished common value in these studies, we have that even though all of them come from different disciplines and topics, in a direct or indirect perspective, they possess the essence of probability and mathematical modeling, and they have an objective directed to solve a particular issue. The works contained in this Research Topic combine all the implicit elements of the general theme, and in the most of them it is mentioned the use of simulation as a primary resource in the handling of studied factors.

Individual works in alphabetical order are:


A Stochastic Phylogenetic Algorithm for Mitochondrial DNA Analysis (Corona-Ruiz et al.):

This research is an exploratory analysis of the mitochondrial DNA and proposes new indices as functions of Shannon entropy, the Chargaff ratio, and fractal dimensions using rescaled-range analysis and DFA; the work suggests to utilize the triplet of indices to construct phylogenetic trees using clustering algorithms. Additionally, it is an initiative to identify the tendencies and correlations in the mutations that produce new species throughout evolutionary history.

On the Aerodynamic Forces on a Baseball, With Applications (Santos et al.):

This article combines experimental results, phenomenological and dimensional analysis; it is classifies as a review, and summarize from recent literature some methods regarding the reproduction and reconstruction of baseball trajectories from aerodynamic forces, and finally discusses their potential applications.

Output-Feedback Control of Virus Spreading in Complex Networks with Quarantine (Alarcn-Ramos et al.):

The research proposes a simple output-feedback control, which stabilizes the extinction state in a virus spreading process over a complex network with quarantine. Additionally, it is provided numerical simulation results to illustrate the functioning of the proposed control scheme for a scale-free network of N = (10)<sup>6</sup> nodes.

Probabilistic Assessment of Investment Options in Honey Value Chains in Lamu County (Wafula et al.):

This work, classified as method brings an analysis to demonstrate a comprehensive approach to decision-making in a project where decision outcomes are uncertain.

The structure of the decision model includes several strategies, of which the Stochastic Impact Evaluation approach (SIE) and Monte Carlo simulation clearly stands out. SIE approach is used as fundamental base for the decision analysis; particularly SIE prioritizes the critical uncertainties considering the risk factors that may compromise the success and the performance of the project. By the other hand, through the use of Monte Carlo simulation, it was possible to obtain the results of project performance over 10 years, to get it possible the model was ran a total of 10,000 times. This propose is substantially useful for the process of decision-making and the probabilistic assessment of investment options, the case of honey value chains in Lamu County, Kenya is analyzed as study case.

Synchronization, Oscillator Death, and Frequency Modulation in a Class of Biologically Inspired Coupled Oscillators (Franci et al.):

This research is focused on the basic mathematical principles that underlie the emergence of synchronous biological rhythms, in particular, the circadian clock. The study analyzes the role that the coupling strength, coupling type, and noise play in the synchronization of a system of coupled, non-linear oscillators. Authors report the following two phenomena found, whose are described as new from a non-linear collective phenomena perspective: In diffusively coupled cells, resting node dynamics imply global asymptotic stability; oscillating node dynamics imply global-synchronization for small coupling, and multistability between oscillator death and global synchronization for large coupling. In stochastic, linearly coupled populations, it is described the dynamical mechanisms through which coupling modulates the frequency of the synchronous oscillation. Finally, the article and its presented results, emphasizes the importance of simple mathematical models in understanding situations where synchronization of multiple oscillating populations appears. Additionally, Authors believe that results obtained may help to shed light on physiological and pathological phenomena involving synchronization of oscillators in important tissues as Parkinson's disease [6, 7], and epilepsy [8, 9].

The Amnesiac Lookback Option: Selectively Monitored Lookback Options and Cryptocurrencies (Chang and Li):

This research suggests the use of its properties to reduce risk exposure in cryptocurrency markets through blockchain enforced smart contracts and correct for informational inefficiencies surrounding prices and volatility. In addition, this work generalizes partial, discretely monitored lookback options that dilute premiums by selecting a subset of specified periods to determine payoff, which authors call amnesiac lookback options. As part of the utilized method, authors price Amnesiac lookbacks with Monte Carlo simulations of Gaussian random walks under equidistant and random periods. The paper concludes that the instrument provides an ideal space for investors to balance their risk, and as a prime candidate to hedge extreme volatility.

## CONCLUDING REMARKS

This collection of articles included in this volume has provided an interesting view of applications associated to fields of probability, statistics, and mathematic modeling, all of them supported by a computational context; additionally it has been determinant the approach of stochasticity and simulation used in most of them as a key element in the utilized methodology. Each article contained, describes a particular object of study, bring a state of art to support the nature of the research, presents a methodology to solve the defined problem and finally provides results and conclusions as the evidence of the efficiency of the entire propose.

We hope the present collection of articles can be attractive for readers from the academic community of mathematics and the applied sciences in general; also, it is our desire that these papers can provide a starting

## REFERENCES


point for future researches in solving of practical and complex problems.

## AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### ACKNOWLEDGMENTS

We want to express our recognition to the Universidad Autónoma de Coahuila and the Universidad Autónoma de Nuevo León for providing all the logistical support to ensure the success of this collaborative work. Likewise, we also express our sincere appreciation to Frontiers and all its team for their kind cooperation and their continuous organizational support, their contribution was key to make possible this project.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Cantú-González and Almaguer. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Output-Feedback Control of Virus Spreading in Complex Networks With Quarantine

Luis A. Alarcón-Ramos <sup>1</sup> , Roberto Bernal Jaquez <sup>2</sup> \* and Alexander Schaum<sup>3</sup>

<sup>1</sup> Posgrado en Ciencias Naturales e Ingeniería, Universidad Autónoma Metropolitana-Cuajimalpa, Mexico City, Mexico, <sup>2</sup> Departamento de Matemáticas Aplicadas y Sistemas, Universidad Autónoma Metropolitana-Cuajimalpa, Mexico City, Mexico, <sup>3</sup> Chair of Automatic Control, Kiel-University, Kiel, Germany

#### Edited by:

Michael Ng, Hong Kong Baptist University, Hong Kong

#### Reviewed by:

Songting Li, New York University, United States Shao-Bo Lin, Wenzhou University, China

#### \*Correspondence:

Roberto Bernal Jaquez rbernal@correo.cua.uam.mx

#### Specialty section:

This article was submitted to Mathematics of Computation and Data Science, a section of the journal Frontiers in Applied Mathematics and Statistics

> Received: 06 April 2018 Accepted: 12 July 2018 Published: 08 August 2018

#### Citation:

Alarcón-Ramos LA, Bernal Jaquez R and Schaum A (2018) Output-Feedback Control of Virus Spreading in Complex Networks With Quarantine. Front. Appl. Math. Stat. 4:34. doi: 10.3389/fams.2018.00034 In this paper the problem of designing an output-feedback control for the stabilization of the extinction steady-state in a virus spreading process over a complex network with quarantine is considered. Sufficient conditions are established for the choice of those nodes for which sensor information is necessary and those which should be controlled using notions from constructive control theory. A simple output-feedback control is proposed which exponentially stabilizes the extinction state. Numerical simulation results are provided to illustrate the functioning of the proposed control scheme for a scale-free network of N = 10<sup>6</sup> nodes.

Keywords: complex networks, virus spreading, feedback control, quarantine, sensor location

#### 1. INTRODUCTION

Studies on the propagation and control of viruses and infectious diseases in human and animal populations [1–8] have gained a great importance in the last years. Understanding and controlling spreading processes is a problem of interdisciplinary nature [9–12]. The use of mathematical and physical inspired models that give account of the dynamics of propagation of infections have gained acceptance since the London cholera epidemics in September 1853 when John Snow divided a map of London into sectors (in a way, nowadays, equiparable to a Voronoi diagram) to calculate who was most likely to use each water pump in the city and in this way, discovered that the 40 Broad Street pump was the main focus of the cholera infection.

Some decades later, McKendrick [13] and Kermack and McKendrick [14] proposed the SIR model that divides the population in three different compartments or groups of individuals: Susceptible, Infective, and Recovered (or Removed) that had the disease and become immune or died. This pioneering work gave birth to the SIS (Susceptible-Infected-Susceptible) model mainly because many diseases do not confer any immunity. Extensions such as the SIQ (Susceptible-Infected-Quarantine) model (see e.g., [15–17]) appeared in order to account for the dynamics of the spreading of infections when quarantine has been implemented as a means to control the massive spreading of the infectious disease. Accordingly, a new class of quarantined individuals is included. These individuals are removed, with some probability from the class of infected individuals.

Although these models were proposed for modeling spreading in populations with no structure, nowadays the SIS and SIQ models have been extended to model the spreading process in a population mapped on a complex network.

In recent years, Markov-chain based models for a Susceptible-Infected-Susceptible (SIS) dynamics over complex networks have been used [4–8, 18] to describe spreading processes in networks. Using these models, it is possible to determine the macroscopic properties of the system as well as the description of the dynamics of individual nodes. One interesting result is that the calculated infection threshold depends on the value of the spectral radius of the adjacency matrix, thus relating the network structure with spreading behavior.

In the present work we will study the spreading of diseases using a Markov-chain based model for a Susceptible-Infected-Quarentine (SIQ) dynamics over a complex network that has been used in a simplified version for the determination of stability tresholds in recent studies by the authors [15, 17]. We will show that, even when the built-in control strategy of quarantine is not able to ensure virus extinction, it is always possible to solve the problem of epidemic spreading extinction. In order to accomplish this task, we determine sufficient conditions for stabilizing the extinction state in a complex network of arbitrary topology. One remarkable point is that, these conditions give a clue to identify the nodes that do not need any control to reach the extinction state and distinguish them from the nodes that need to be controlled. Inspired on the ideas of control theory, the set of nodes that do not need to be controlled (in order to reach the extinction state) will be associated with the zero dynamics of the system [19, 20] . At the same time, the set of nodes to be monitored and controlled will be identified and a decentralized feedback control will be applied in order to stabilize the extinction state. Accordingly, it is proven that the extinction state is an exponentially stable fixed point for the zero dynamics.

We have performed numerical simulations using a scalefree complex network with 1 million nodes constructed as proposed by Barabasi [21] and found a complete agreement of the numerical results with our theoretical findings.

The paper is organized as follows: section Methods presents the problem statement for the SIQ model mapped on a complex network of arbitrary topology and the control problem we want to solve together with some definitions. In section Results, sufficient conditions for the stabilization of the extinction state are derived. Using these conditions we establish a selection criterion that allows to identify the set of nodes that need to be controlled in order to reach the extinction state. Afterwards, we design a simple stabilizing output-feedback control and present our simulation results. In section Discussion and Conclusions we summarize our main results and present our main conclusions.

#### 2. METHODS

As pointed out in the introduction, the SIQ model belongs to a class of models that are used to capture possible human interactions employed to impede the spreading of an infection process and it is essentially an extension of the SIS model in which the model and control strategies are co-developed to yield a kind of closed-loop control model [22]. In order to reveal the dynamics and essential mechanisms of this model, we proceed to (i) formulate the SIQ model in a complex network (ii) introduce the specific control problem.

## 2.1. The SIQ Model in a Complex Network

Consider a network of N nodes described by an undirected graph G(V, E) of any topology. Let V = {v1, v2, . . . , vN} being the set of nodes and E = {ei,j} the set of connecting edges. The adjacency matrix associated to G(V, E) is given by A = {aij}, where aij = 1 if ei,<sup>j</sup> ∈ E and zero otherwise. The set of neighbor nodes of a node v<sup>i</sup> ∈ V is defined as

$$V\_i = \{ \nu\_j \in V \mid a\_{ij} = a\_{ji} = 1 \} \subset V,\tag{1}$$

and the number of neighbors or degree of a node v<sup>i</sup> is given by N<sup>i</sup> = |V<sup>i</sup> |.

As formulated in Bernal Jaquez et al. [15], the underlying process for every node is depicted in **Figure 1**, as a discrete time Markov process. A node v<sup>i</sup> can be in state I (Infected) with probability pi(t), in state Q (Quarantine) with probability qi(t), or in state S (Susceptible) with probability si(t) = 1 − pi(t) − qi(t). At each time step, the probability functions pi(t),si(t) are updated due to the fact that every node can transit from state S to state I with probability 1 − ζ<sup>i</sup> or from state I to state Q with probability τ<sup>i</sup> and from state Q to state S with probability µ<sup>i</sup> because nodes are interacting.

According to the transition diagram shown in **Figure 1**, we have the following dynamical system:

$$\begin{aligned} p\_i(t+1) &= (1 - \tau\_i)p\_i(t) + (1 - \xi\_i(P(t)))s\_i(t), \\ q\_i(t+1) &= (1 - \mu\_i)q\_i(t) + \tau\_i p\_i(t), \\ s\_i(t+1) &= \xi\_i(P(t))s\_i(t) + \mu\_i q\_i(t). \\ p\_i(t) + q\_i(t) + s\_i(t) &= 1, \quad \nu\_i \in V, \end{aligned} \tag{2}$$

where, for every node v<sup>i</sup> ∈ V, τ<sup>i</sup> is the internment probability associated to quarantine, µ<sup>i</sup> is the recovery probability and

ζi(P(t)) is the probability of a node v<sup>i</sup> of not being infected by any neighbor at time t, which is given by

$$\langle \zeta\_i(P(t)) \rangle = \prod\_{j=1}^{N} (1 - \beta\_i r\_{ij} a\_{ij} p\_j(t)). \tag{3}$$

In the above equation, we have defined the vector P(t) = [p1(t), p2(t), ..., pN(t)]<sup>T</sup> ∈ [0, 1]<sup>N</sup> and β<sup>i</sup> is the probability of infection during a single contact, rij is the probability that the node v<sup>i</sup> performs at least one contact intent with its neighbor v<sup>j</sup> ∈ V<sup>i</sup> . The probability rij is also known as the connection probability and depends on the number of contact intents (or interaction rate).

Note that in system (2)

$$0 \le \mu\_i, \ r\_i, r\_{ij}, \beta\_i, \zeta\_i(P(t)), p\_i(t), q\_i(t), s\_i(t) \le 1$$

From the model (2) we point out that


We consider that each node has a manipulable variable ui(t), which is amenable for control. For the system 2, this manipulable variable can be β<sup>i</sup> or rij, i.e., we consider that, for each node v<sup>i</sup> , it is possible to improve its health or avoid to perform several contact attempts with its neighbor nodes, the control will adapt one of these parameters, that will be selected according to our mathematical analysis.

#### 2.2. Control Problem

Quarantine is an heuristic control mechanism or strategy that intents to reduce a virus propagation in a population. As stated in the last subsection, these kind of models probably are not sufficient to ensure the system to reach the extinction state. This is due to the fact that no adaptation of the network parameters is implemented in dependence of the actual system state. Such a mechanism is proposed in the sequel including decisions on (i) the subset of nodes V<sup>M</sup> ⊂ V which should be monitored, i.e., whose actual state must be known at each time instant, (ii) the subset of nodes V<sup>C</sup> ⊂ V which have to be controlled, i.e., whose interaction parameters (u<sup>i</sup> = β<sup>i</sup> or u<sup>i</sup> = rij) should be subject to on-line adaptation, and (iii) the specific control law ui(t) = ϕi(P(t)) which should be used for this adaptation.

The chosen approach follows the constructive (i.e., passivity-based) control idea, and consist in two steps:


$$p\_i(t) \le p\_{io} \nu^t, \quad p\_i(t=0) = p\_{io} \dots$$

#### 3. RESULTS

In this section the main results are presented. In particular sufficient conditions for the stabilization of the extinction state are derived including (i) a selection criterion for the nodes to be monitored in terms of the connectivity parameters, the infection probability and the graph topology and (ii) the design of a simple stabilizing output-feedback control scheme. Simulation results illustrate the functioning of the proposed control scheme for a scale-free network with N = 10<sup>6</sup> nodes.

#### 3.1. Selection of Monitored and Controlled Nodes

Taking into account that pi(t) + qi(t) + si(t) = 1 for each v<sup>i</sup> ∈ V and to be consistent with the control idea, we consider that ζi(P) = ζi(P, U), where the vector U(t) = [uc1(t), uc2(t), . . . , ucK(t)]<sup>T</sup> , vci ∈ V<sup>C</sup> represent the manipulable set of parameters associated with β<sup>i</sup> or rij. So, we can rewrite system (2) as follows

$$\begin{aligned} p\_i(t+1) &= (1 - \tau\_i)p\_i(t) + (1 - \xi\_i(P(t), U(t)))(1 - p\_i(t) - q\_i(t)), \\ q\_i(t+1) &= (1 - \mu\_i)q\_i(t) + \tau\_i p\_i(t), \\ \xi\_i(P, U) &= \prod\_{j=1}^{N} (1 - \beta\_i r\_{ij} a\_{ij} p\_j(t)). \end{aligned} \tag{4}$$

The fixed points associated with the dynamics (4) for some constant U ∗ (i.e., β<sup>i</sup> and rij are set to some constant value) can be determined by substituting the relations pi(t + 1) = pi(t) = p ∗ i and qi(t + 1) = qi(t) = q ∗ i . After some algebra it follows that

$$\begin{split} p\_i^\* &= \frac{1 - \zeta\_i(P^\*, U^\*)}{\tau\_i + \left(1 + \frac{\tau\_i}{\mu\_i}\right) \left(1 - \zeta\_i(P^\*, U^\*)\right)}, \\ q\_i^\* &= \frac{\tau\_i}{\mu\_i} \frac{1 - \zeta\_i(P^\*, U^\*)}{\tau\_i + \left(1 + \frac{\tau\_i}{\mu\_i}\right) \left(1 - \zeta\_i(P^\*, U^\*)\right)}, \\ s\_i^\* &= \frac{\tau\_i}{\tau\_i + \left(1 + \frac{\tau\_i}{\mu\_i}\right) \left(1 - \zeta\_i(P^\*, U^\*)\right)}. \end{split} \tag{5}$$

Note that the extinction state (p ∗ <sup>i</sup> = 0, q ∗ <sup>i</sup> = 0,s ∗ <sup>i</sup> = 1) in (2), for v<sup>i</sup> ∈ V is a fixed point when ζi(0, U ∗ ) = 1. However, up to this point, it is not clear if the extinction state or any other fixed point given by (5) are stable. The extinction state means that no viruses are propagated over the network and, as we will see in the following Lemma, knowledge of the conditions under which this state is reached, will give us a clue of how viruses are propagated.

The approach followed subsequently exploits ideas from Wang et al. [4], and Bernal Jaquez et al. [15] by establishing a linear bounding dynamics for (4) that has the origin (P, Q) = (0, 0) as exponentially stable fixed point.

Lemma 1. Consider the dynamics (4) on a complex network with graph G(V, E) and adjacency matrix **A**. The extinction state (P, Q) = (0, 0) is globally exponentially stable if the constant vector U ∗ (i.e., for some constant values of β<sup>i</sup> and rij) is such that

$$
\sigma(\mathbf{H}) < 1, \quad U^\* = [\mu\_1^\*, \dots, \mu\_N^\*], \tag{6}
$$

where σ(·) is the spectral radius of the matrix **H** defined as

$$\mathbf{H} = \begin{bmatrix} \mathbf{I} - \mathbf{T} + \mathbf{BR} & \mathbf{0} \\ \mathbf{T} & \mathbf{I} - \mathbf{M} \end{bmatrix},\tag{7}$$

where

$$\mathbf{T} = \operatorname{diag}(\mathbf{r}\_i), \quad \mathbf{B} = \operatorname{diag}(\beta\_i), \quad \mathbf{R} = [r\_{ij} a\_{ij}], \quad \mathbf{M} = \operatorname{diag}(\mu\_i). \tag{8}$$

and **I** being the identity matrix.

Proof: As it is proved in Wang et al. [4], and Bernal Jaquez et al. [15], ζi(P, U ∗ ) can be bounded as follows

$$1 - \zeta\_i(P(t), U^\*) \le \sum\_{j=1}^N \beta\_i r\_{ij} a\_{ij} p\_j(t). \tag{9}$$

Substituting this bound into the first equation in (4) and after some algebra one obtains

$$\begin{aligned} p\_i(t+1) &\le (1 - \tau\_i)p\_i(t) + (1 - p\_i(t) - q\_i(t)) \sum\_{j=1}^N \beta\_i r\_{ij} a\_{ij} p\_j(t), \\ &\le (1 - \tau\_i)p\_i(t) + \sum\_{j=1}^N \beta\_i r\_{ij} a\_{ij} p\_j(t). \end{aligned}$$

This can be written in matrix form as

$$
\begin{bmatrix} P(t+1) \\ Q(t+1) \end{bmatrix} \le \begin{bmatrix} \mathbf{I} - \mathbf{T} + \mathbf{BR} & \mathbf{0} \\ \mathbf{T} & \mathbf{I} - \mathbf{M} \end{bmatrix} \begin{bmatrix} P(t) \\ Q(t) \end{bmatrix},
$$

where P(t) = [p1(t), . . . , pN(t)]<sup>T</sup> and Q(t) = [q1(t), . . . , qN(t)]<sup>T</sup> , and the inequality being interpreted as element-wise. Therefore, the solutions of pi(t) are bounded by the linear dynamics

$$\begin{aligned} \mathbf{x}\_{i}(t+1) &= (1 - \mathbf{r}\_{i})\mathbf{x}\_{i}(t) + \sum\_{j=1}^{N} \beta\_{i} r\_{ij} a\_{ij} \mathbf{x}\_{j}(t), \\ \mathbf{w}\_{i}(t+1) &= \mathbf{r}\_{i} \mathbf{x}\_{i}(t) + (1 - \mu\_{i})\mathbf{w}\_{i}(t), \quad \mathbf{v}\_{i} \in V \end{aligned} \tag{10}$$

i,e., forall t ≥ 0 it holds that

$$0 \le p\_i(t) \le \varkappa\_i(t), \quad q\_i(t) = \varkappa\_i(t) \tag{11}$$

if p(0) = x(0), q(0) = w(0). The linear dynamics can be expressed in matrix form as

$$
\begin{bmatrix} X(t+1) \\ W(t+1) \end{bmatrix} = \begin{bmatrix} \mathbf{I} - \mathbf{T} + \mathbf{BR} & \mathbf{0} \\ \mathbf{T} & \mathbf{I} - \mathbf{M} \end{bmatrix} \begin{bmatrix} X(t) \\ W(t) \end{bmatrix} =: \mathbf{H} \begin{bmatrix} X(t) \\ W(t) \end{bmatrix}, \tag{12}
$$

where X(t) = [x1(t), . . . , xN(t)]<sup>T</sup> and W(t) = [w1(t), . . . ,wN(t)]<sup>T</sup> . It holds that (X T , W<sup>T</sup> ) <sup>T</sup> = (0<sup>T</sup> , 0<sup>T</sup> ) T is exponentially stable if and only if the eigenvalues of the associated matrix **H** are contained in the open unit circle C<sup>1</sup> = {λ ∈ C | |λ| < 1}. Taking into account (11) it follows that a necessary and sufficient condition for global exponential stability of (P T , Q T ) <sup>T</sup> = 0 is given by (6).

The expontial stability condition (6) is very general and does not provide any idea on how to select nodes to be controlled or monitored in order to reach the extinction state. However, using this result as a point of departure, we can get insight into the condition that every node v<sup>i</sup> ∈ V has to fulfill in order to ensure that P(t) converges to P <sup>∗</sup> = 0 as shown in the following Lemma.

Lemma 2. For a constant U<sup>∗</sup> (i.e., for some constant value for β<sup>i</sup> and rij), the state vector (P T (t), Q T (t))<sup>T</sup> globally exponentially converges to the extinction state (P T , Q T ) <sup>T</sup> = 0 if for every node v<sup>i</sup> ∈ V it holds that

$$
\beta\_i < \frac{r\_i}{\sum\_{j=1}^N r\_{ij} a\_{ij}}.\tag{13}
$$

Proof: In virtue of Lemma 1 it is sufficient to show that if the condition (13) holds, the matrix **H** defined in (6) has only eigenvalues within the open unit circle C<sup>1</sup> ⊂ C, or equivalently, that its spectral radius σ(**H**) < 1. Note that due to the blockdiagonal structure of the matrix **H** its eigenvalues are given by the eigenvalues of the two matrices on its diagonal, i.e., in terms of the matrix spectra

$$\mathcal{S}(\mathbf{H}) = \mathcal{S}(\mathbf{I} - \mathbf{T} + \mathbf{BR}) \cap \mathcal{S}(\mathbf{I} - \mathbf{M}),$$

where S(**A**) denotes the spectrum of a matrix **A**, i.e., the union of all its eigenvalues. The eigenvalues of the matrix **I** − **M** are contained in C1, given that **M** is diagonal with entries 0 < µ<sup>i</sup> < 1. Thus ,it remains to show that S(**I** − **T** + **BR**) ⊂ C1. This can be analyzed by applying Gerschgorin's Theorem [23], which provides an upper-bound estimate for the spectral radius of a given matrix. In the following, let λ represent an arbitrary eigenvalue of the matrix **I** − **T** + **BR**. The application of this theorem to the matrix **I** − **T** + **BR** provides the following inequality

$$|\lambda| \le |1 - \mathfrak{r}\_i + \beta\_i \sum r\_{ij} a\_{ij}| = 1 - \mathfrak{r}\_i + \beta\_i \sum r\_{ij} a\_{ij} \tag{14a}$$

i = 1, . . . , N. Thus, |λ| < 1 is satisfied if

$$\sum\_{j=1}^{N} \beta\_i r\_{ij} a\_{ij} < \mathfrak{r}\_i, \quad \forall \mathfrak{v}\_i \in V$$

or equivalently if (13) holds true. This complete the proof.

Note that, as stated in the proof the convergence to zero of X (and thus P) is independent of the behavior of W (or accordingly Q). This resides in the fact that the dynamics represent a cascade structure. Accordingly, for the purpose of stabilizing the extinction state it is sufficient to ensure that less nodes get infected than pass into quarantine. This intuitively clear condition is exactly what is formally stated in Lemma 2.

The condition (13) gives a criterion on how to choose the nodes to be monitored. If condition (13) does not hold for some set of nodes, then it is appropriate to consider this collection as the set of nodes to be monitored VM. Additionally, according to (13), the set of controlling nodes can be set as V<sup>C</sup> = VM. That is, we can consider that all controlled nodes are monitored. In this way, we can select β<sup>i</sup> (or rij) as the parameter amenable for control purposes for those nodes that do not satisfy (13).

Further note that in case that the network parameters are homogenous, i.e., all nodes have the same parameter values µ, β, τ ,rij = r, the selection criterion for a node to be controlled is directly related to its degree N<sup>i</sup> = P<sup>N</sup> i=1 aij (i.e., here simply the number of neighboring nodes) and can be expressed as

FIGURE 3 | Behavior of ρ(t) for several initial conditions without control.

v<sup>i</sup> ∈ V<sup>C</sup> ∀ i : N<sup>i</sup> > τ rβ . Nevertheless, when the parameter distribution is nonhomogenous, it is possible that nodes with high degree do not have to be controlled and nodes with small degree have to be, given the particular constellation between τ<sup>i</sup> , β<sup>i</sup> and rij according to condition (13). This is illustrated in **Figure 2** for a scale-free network with N = 100 nodes and normally distributed parameters.

Note that in the notion of constructive control theory (see e.g., [19, 24]), the dynamics of the nodes that are not controlled establishes the zero dynamics (i.e., the dynamics resulting from the restriction that p<sup>i</sup> = 0 for all i such that v<sup>i</sup> ∈ Vc) and by construction this dynamics is asymptotically stable. Therefore, the zero dynamics correspond to a spreading process over a reduced network, from which the monitored (and thus controlled) nodes have been withdrawn, given that for pi(t) ≡ 0 the node v<sup>i</sup> does no longer interact with its neighbors.

#### 3.2. Feedback Control Design

The question addressed in this section is how to design the feedback control for the nodes v<sup>j</sup> ∈ V<sup>C</sup> so that limt→∞ |pj(t) − p ∗ j | = 0. Up to this point, we have not considered the dependency of ζ (P, U) on U, given that U was considered as a set of constant parameters. This dependency will permit to explicitly determine a control law ui(t) = βi(t) that steers the nodes v<sup>j</sup> ∈ V<sup>C</sup> to their desired values p ∗ <sup>j</sup> = 0. With this idea, the following theorem is established

Theorem 1. Consider the dynamics given by (4) and consider the nodes that do not satisfy (13) as the set V<sup>C</sup> of nodes to be controlled. Let V<sup>M</sup> = VC, that is, all controlled nodes are monitored. If the controls ui(t) satisfy

$$0 \le \mu\_i(t) < \bar{\beta}\_i = \frac{\tau\_i}{\sum\_{j=1}^N r\_{ij} a\_{ij}},\tag{15}$$

then (P T , Q T ) <sup>T</sup> = 0 is exponentially stable.

Proof: The Theorem can be easily proven as follows. Let β¯ <sup>i</sup> be given as in (15). It follows that

$$\begin{aligned} 0 \le \mu\_i(t) < \bar{\beta}\_i, \\ \Leftrightarrow \quad 0 \le \mu\_i(t)r\_{ij}a\_{ij} < \bar{\beta}\_i r\_{ij}a\_{ij}, \end{aligned}$$

$$\Leftrightarrow \quad 0 \le \sum\_{j=1}^{N} \mu\_i(t) r\_{ij} a\_{ij} < \sum\_{j=1}^{N} \bar{\beta}\_i r\_{ij} a\_{ij}.$$

Given that for β<sup>i</sup> ≤ β¯ i inequality (13) is satisfied, one obtains

$$0 \le \sum\_{j=1}^N u\_i(t)r\_{ij}a\_{ij} < \sum\_{j=1}^N \bar{\beta}\_i r\_{ij}a\_{ij} < \mathfrak{r}\_i,$$

or equivalently

$$\sum\_{j=1}^{N} u\_i(t) r\_{ij} a\_{ij} < \tau\_i$$

implying that inequality (13) is satisfied. In consequence, the exponential stability of (P T , Q T ) <sup>T</sup> = 0 follows.

Note that following a similar reasoning, a control law can be established for the case that ui(t) = rij(t). Nevertheless, this approach will not be explicitely elaborated at this place.

From Theorem 1, we have that in order to stabilize the extinction state for the dynamics (4) it is sufficient to design feedback controls ui(t) for all nodes v<sup>i</sup> ∈ VC, which take values below the upper-bound β¯ <sup>i</sup> defined in (15). For example, a simple linear feedback control given by

$$u\_i(t) = \beta\_i(t) = \gamma \,\vec{\beta}\_i (1 - p\_i(t)), \quad \gamma \in (0, 1) \tag{16}$$

does satisfy this condition. Note that the advantage of using a time varying control law βi(t) consists in actively adapting the infection probability on the actual needs, i.e., the actual network state. In comparison with imposing a constant value for β<sup>i</sup> this possibly enables to optimize the control effort over time, because according to (16) β<sup>i</sup> → γ β¯ <sup>i</sup> with p<sup>i</sup> → 0. The performance using this simple output-feedback control is illustrated in the subsequent section.

#### 3.3. Simulations

In order to verify our results, we perform several simulations of the dynamical system (4), for different initial conditions, with the following considerations


• To show and corroborate our results, we calculated the average probability ρ(t) given by

$$\rho(t) = \frac{1}{N} \sum\_{i=1}^{N} p\_i(t). \tag{17}$$

• The sets V<sup>C</sup> and V<sup>M</sup> are chosen according to Lemma 2 as the sets of nodes that do not satisfy (13). For the case considered here 424569 nodes out of N = 10<sup>6</sup> (i.e., about 42.5%) have to monitored and controlled.

Simulation studies have been performed considering the following scenarios:


#### I. Absence of Control

In the absence of control the network reaches an endemic attractor state where about a 20% of the nodes are probably infected after about 20 time units as can be seen in **Figure 3**.

#### II. Zero Dynamics

To corroborate that the zero dynamics are exponentially stable, we set Y(t) = 0 for those nodes that do not satisfy (13). This is achieved by setting β<sup>i</sup> = 0 for v<sup>i</sup> ∈ VM.

**Figure 4** shows the behavior of the nodes associated with the zero dynamics. Note that the extinction state is reached after 18 time steps.

#### III. Linear Feedback Control

To show the performance of the proposed simple outputfeedback control scheme (16), it has been implemented for all monitored and controlled nodes with the gain γ = 0.9 and the upper bound β¯ <sup>i</sup> calculated according to (15). This yields the output-feedback controller.

$$\beta\_i(t) = \frac{\boldsymbol{\gamma} \,\tau\_i (1 - p\_i(t))}{\sum\_{j=1}^{N} r\_j a\_{ij}}, \quad \text{for} \quad \boldsymbol{\gamma} = 0.9 \quad \text{and} \quad \forall \boldsymbol{\nu}\_i \in V\_C = V\_M. \tag{18}$$

#### REFERENCES


Note that the control (18), depends on the state of the node i given by pi(t), and the properties of its neighbors given by β¯ i . **Figure 5** shows the result of the simulation with the applied control (18). As predicted in Theorem 1 the extinction state is a close-loop attractor, and is reached in about 40 time steps.

#### 4. DISCUSSIONS AND CONCLUSIONS

The problem of deciding which nodes in a complex network with quarantine should be controlled and how to control them in order to achieve that in the closed-loop system the extinction state becomes a (global) attractor has been studied. Sufficient conditions for virus extinction have been derived using a constructive control approach by suitably identifying the zero dynamics according to a threshold condition for the stability of the extinction state. The associated node selection criterion does depend on the transmission probability between the nodes and their neighbors, the degree of each node and its probability to pass into quarantine. It has been shown that in spite of the strongly nonlinear dynamics of the spreading process the extinction state can be efficiently stabilized using a simple linear bounded output-feedback control if the nodes to be controlled are selected according to the proposed scheme. The performance and behavior of the spreading process without and with control has been illustrated for a scale-free network with N = 10<sup>6</sup> nodes.

#### AUTHOR CONTRIBUTIONS

All authors designed and did the research, in particular the stability analysis. LA-R and AS conceived and designed the control and the numerical experiments. LA-R performed the numerical experiments. All authors analyzed the data. LA-R and RB wrote the paper. All authors reviewed the manuscript. All authors have read and approved the final manuscript.

#### FUNDING

This work was partially sponsored by Rector office of UAM Cuajimalpa through Programa de Apoyo a Proyectos Interdisciplinarios.


and non-homogeneous transition rates. J Phys Conf Ser. (2014) 490:012011. doi: 10.1088/1742-6596/490/1/012011


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Alarcón-Ramos, Bernal Jaquez and Schaum. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Synchronization, Oscillator Death, and Frequency Modulation in a Class of Biologically Inspired Coupled Oscillators

Alessio Franci <sup>1</sup> \*, Marco Arieli Herrera-Valdez <sup>1</sup> , Miguel Lara-Aparicio<sup>1</sup> and Pablo Padilla-Longoria<sup>2</sup> \*

<sup>1</sup> Departamento de Matemáticas, Facultad de Ciencias, Universidad Nacional Autónoma de México, Mexico City, Mexico, <sup>2</sup> Departamento de Matemáticas y Mecánica, Instituto de Investigación en Matemáticas Aplicadas y Sistemas, Universidad Nacional Autónoma de Mexico, Mexico City, Mexico

#### Edited by:

José Roberto Cantú-González, Universidad Autónoma de Coahuila, Mexico

#### Reviewed by:

Alessandro Giuliani, Istituto Superiore di Sanità (ISS), Italy Alexey Goltsov, Abertay University, United Kingdom

#### \*Correspondence:

Alessio Franci afranci@ciencias.unam.mx Pablo Padilla-Longoria pablo@mym.iimas.unam.mx

#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Applied Mathematics and Statistics

> Received: 07 May 2018 Accepted: 12 October 2018 Published: 08 November 2018

#### Citation:

Franci A, Herrera-Valdez MA, Lara-Aparicio M and Padilla-Longoria P (2018) Synchronization, Oscillator Death, and Frequency Modulation in a Class of Biologically Inspired Coupled Oscillators. Front. Appl. Math. Stat. 4:51. doi: 10.3389/fams.2018.00051 The general purpose of this paper is to build up on our understanding of the basic mathematical principles that underlie the emergence of synchronous biological rhythms, in particular, the circadian clock. To do so, we study the role that the coupling strength, coupling type, and noise play in the synchronization of a system of coupled, non-linear oscillators. First, we study a deterministic model based on Van der Pol coupled oscillators, modeling a population of diffusively coupled cells, to find regions in the parameter space for which synchronous oscillations emerge and to provide conditions under which diffusive coupling kills the synchronous oscillation. Second, we study how noise and coupling interact and lead to synchronous oscillations in linearly coupled oscillators, modeling the interaction between various pacemaker populations, each having an endogenous circadian clock. To do so, we use the Fokker-Planck equation associated to the system. We show how coupling can tune the frequency of the emergent synchronous oscillation, which provides a general mechanism to make fast (ultradian) pacemakers slow (circadian) and synchronous via coupling. The basic mechanisms behind the generation of oscillations and the emergence of synchrony that we describe here can be used to guide further studies of coupled oscillations in biophysical non-linear models.

Keywords: synchronization, oscillator death, circadian rhythm, coupled non-linear oscillators, Fokker-Planck

## 1. INTRODUCTION

The study of circadian rhythms has been a subject of great interest for a long time. The majority of the first studies were mainly based on observations in plants [1–4]. The study of circadian rhythms from a mathematical perspective reached a milestone with the work of Kalmus and Wigglesworth, a biologist and a mathematician, respectively, who associated of circadian rhythms to the existence of a limit cycle, using a hydraulic system as analogy. Kalmus and Wigglesworth presented their work entitled "Shock excited system as models for biological rhythms" along with several mathematical models of circadian rhythms at a Symposium on Biological Rhythms carried out at Cold Spring Harbor in the United States of America in 1960 [5–7]. Lots of other important works on circadian rhythms were presented in this symposium, but the work of Kalmus and Wigglesworth was key in establishing a better mathematical formalism for the study of circadian rhythms. Although many researchers followed the theoretical path proposed by Kalmus and Wigglesworth, the mathematical study of circadian rhythms was finally established by Arthur Winfree (biologist and mathematician), who introduced topology for the description of several aspects of circadian rhythms. An excellent summary of many of the earlier works can be found in Arthur Winfree's master book entitled "The Geometry of Biological Time" [8]. The number of studies about biological rhythms at large has increased greatly in the last two decades, in part due to new technological advances. Particularly related to circadian rhythms, there is now evidence that there is rhythmic patterns of activity at the molecular [9, 10], cellular [11, 12], tissue [13, 14], and systems levels [15–17], and that circadian regulation is involved in jointly regulating activity in all those different levels of biological organization [18–20], and also, taking into account interactions with the environment [21] and perturbations induced by behavior [22, 23].

Mathematical modeling and experimental characterizations of different properties of circadian rhythms have been combined to produce explanations and hypotheses about the rhythmicity in biological phenomena [24–28]. Of particular interest, the ontogeny of circadian rhythm in the crayfish has been studied by Lara-Aparicio et al. [29] combining theoretical and experimental perspectives, building phenomenological mathematical models that capture a series of experimental results involving the synchronization of electro-retinogram activity in crayfish [30– 32]. These models couple two van der Pol oscillators [33] represented by state vectors, and include parameters representing the frequency of the oscillators, the radii of the limit cycles, and the first coordinate of the center of the limit cycle. The system is setup such that one oscillator is driving the behavior of the other oscillator, but not vice-versa. One of the main findings with this model is that the driving oscillator induces an Andronov-Hopf bifurcation in the driven oscillator and regulates its frequency.

The model by Aparicio et al. simulates, explains, and has suggested new biological experiments, it is simple enough to allow analytical approaches, and it has provided useful insights about questions referring to the ontogeny of the circadian rhythm in crayfish from the childhood to adult stages. For instance, a hypothesis about the existence of a hormone, which was experimentally detected, was generated from the model. The model also allowed Lara-Aparicio et al. to study synchronization of circadian rhythms with external signals like day and night light cycle [34]. By studying basic principles underlying the generation of oscillations in coupled non-linear systems, these researchers were able to conjecture that circadian rhythms can result from coupling systems of cells, each one oscillating with an ultradian oscillation [35–37]. Synchronization among cells emerges naturally as a motivating theme that has been studied through systems of non-linear Equations [38] representing n oscillators with the classical van der Pol non-linear damping for the terms responsible for the oscillatory dynamics.

In the present paper, we extend the work in Lara-Aparicio et al. [29] and Barriga-Montoya et al. [38] by analyzing two qualitative mathematical models, each capturing a different level of organization in the ontogeny of circadian rhythms. Inspired by gap junction coupling between neurons, or similarly, by chemical coupling in self oscillating networks, we study the bifurcation structure in a deterministic model assuming that the coupling between the oscillators is diffusive. The resulting dynamics resemble neuronal activity at the cellular level. Then, using graph theoretical methods and center manifold theory, we show that synchronous oscillations appear via a Hopf bifurcation in a population of pacemaker oscillators. In this case, the bifurcation parameter is thought of as a representative of the developmental stage of the neurons. We further explore the phenomenon of oscillator death: although the single neurons are endogenously oscillating for sufficiently large values of the bifurcation parameter, the population oscillation dies for sufficiently large coupling, which suggests that the weak coupling hypothesis must be satisfied for robust synchronous oscillations to occur. In the case of all-to-all coupling, we provide necessary conditions for oscillator death to occur and leave the derivation of sufficient conditions to a future report.

Frequency modulation is also an important phenomenon that can be studied with these models. In consideration of the results from the first analysis, we shift our focus to the frequency modulations that emerge as a result of the interconnection of various circadian pacemaker populations. In doing so, we estimate the synchronization frequency as a function of the intrinsic frequencies of the oscillators, their coupling strength, and the topology of the network. To do so, we construct a second model that can be thought of as a stochastic, larger-scale extension of the first model we present. In this case, each population is assumed to be an endogenous oscillator and the coupling is assumed to be linear. Using linear stochastic analysis and under the assumption that the population oscillations are synchronized, we derive a scalar Fokker-Plank equation. The model captures an important feature of circadian rhythm ontogeny: the emergence of low frequency (circadian) oscillations from coupled high frequency (ultradian) oscillators [35–38]. Future work will aim at deriving conditions on the intrinsic frequencies, the coupling strength, and the network topology, that ensure synchronization of the population oscillations.

It is worth noticing that, although the motivation for the present paper is centered around circadian rhythms, the model captures more general phenomena. Our findings include that in diffusively coupled cells resting node dynamics imply global asymptotic stability, oscillating node dynamics imply global synchronization for small coupling, and multistability between oscillator death and global synchronization for large coupling. In stochastic, linearly coupled populations, we describe the dynamical mechanisms through which coupling modulates the frequency of the synchronous oscillation. To the author's knowledge, both phenomena are new from a non-linear collective phenomena perspective. Among other reasons, it is surprising that passive coupling like diffusive coupling could kill an oscillation and create multistability. Similarly, we are not aware of any work providing mechanistic explanations on how coupling can tune a global oscillation frequency.

The paper is organized as follows. In section 2, we present and analyze the model of diffusively coupled oscillators. Two theorems are proved about global asymptotic stability for resting cells and global synchronization for oscillating cells and weak coupling. We then derive sufficient conditions for diffusive coupling-induced oscillator death. In section 3 we present and analyze the model of linearly coupled oscillators. In particular, we derive an explicit formula for the emergent synchronization frequency as a function of the coupling topology and oscillator natural frequencies. Finally we discuss the presented results in section 4.

## 2. GLOBAL SYNCHRONIZATION AND OSCILLATOR DEATH IN DIFFUSIVELY COUPLED OSCILLATORS

We regard a network of N coupled oscillators as a directed graph G with N vertices, with a network topology codified by a matrix A = [aij] N i,j=1 , where aij ≥ 0 represent connection weights. If oscillator i receives signals from oscillator j, then the graphical representation of G has an arrow from j to i, and aij > 0. If aij > 0, the signal received by oscillator i from oscillator j is aij(x<sup>j</sup> − xi). Assume that the dynamics for each oscillator satisfies the following coupled oscillator dynamics

$$\dot{\boldsymbol{x}}\_{i} = \boldsymbol{\gamma}\_{i} + \mu \sum\_{j=1}^{N} a\_{ij} (\boldsymbol{x}\_{j} - \boldsymbol{x}\_{i}) \tag{1a}$$

$$
\dot{\nu}\_i = \left(\lambda - \mathbf{x}\_i^2 - \frac{\mathbf{y}\_i^2}{\alpha^2}\right)\nu\_i - \alpha\_i^2 \mathbf{x}\_i \tag{1b}
$$

where µ is the global coupling strength, which uniformly scales the coupling weights aij.

In the absence of coupling the oscillator dynamics reduces to the modified van der Pol equation

$$
\ddot{\boldsymbol{x}}\_i = \left(\lambda - \boldsymbol{x}\_i^2 - \frac{\boldsymbol{\nu}\_i^2}{\boldsymbol{\alpha}^2}\right) \dot{\boldsymbol{x}}\_i - \boldsymbol{\alpha}\_i^2 \boldsymbol{\kappa}\_i
$$

These equations define a simple dynamical system that can transition between global asymptotic stability and almost global convergence to a hyperbolic limit cycle through variations of the control parameter λ ∈ R, providing a simple model for various biological systems that exhibit the same transition, in particular neurons and molecular oscillators. For λ < 0 the nonlinear dissipation coefficient − λ − x <sup>2</sup> − y 2 ω 2 i is always positive, which leads to damped oscillations. For λ > 0 the dissipation coefficient becomes negative close to the origin, which leads to sustained oscillations through a Hopf bifurcation. A generic trajectory belonging to the family of periodic orbits born at the Hopf bifurcation has the form

$$
\sqrt{\lambda} \langle \cos(\omega\_i t + \theta\_0), \sin(\omega\_i t + \theta\_0) \rangle,\tag{2}
$$

where the θ<sup>0</sup> is the initial phase.

In the presence of coupling, Equations (1) represent a generic network of diffusively coupled non-linear oscillators. As mentioned earlier, the diffusive form of the coupling can be thought of as gap junction coupling in a neuronal population, or diffusive coupling between non-linear molecular oscillators. In both interpretations, the graph topology is necessarily undirected, that is, aij = aji. However, the mathematical results presented in this section hold under more general assumptions and we present them in the general form.

#### 2.1. Global Synchronization

We start by recalling some basic graph-theoretical definitions and facts. The graph G is said to be strongly connected if for each pair of nodes in G, there exists a directed path between them. G is balanced if P<sup>N</sup> j=1 aij = P<sup>N</sup> j=1 aji for all i.

The Laplacian matrix L = - Lij : i, j = 1, ..., n for the graph G is such that Lij = −aij if i 6= j, and Lii = P<sup>N</sup> j=1 aij. Note that the vector of ones is always a right null eigenvector of L and zero is always an eigenvalue of L (L1<sup>N</sup> = 0). It can be shown that a graph is strongly connected if and only if zero is a simple eigenvalue of the Laplacian matrix [39]. Obviously, symmetric graphs (i.e., satisfying aij = aji) are balanced, but the converse is not true. Consider, for instance, a directed ring.

The global behavior of the system (1) for λ < 0 and µ ≥ 0, for a network with connectivity represented by a generic balanced graph is characterized by the next theorem (**Figure 1**).

**Theorem 2.1.** Assume that the graph G is balanced and that ω<sup>i</sup> = ω<sup>j</sup> = ω for all i, j = 1, . . . , N. If λ < 0, then the origin is globally asymptotically stable and locally exponentially stable for all µ ≥ 0.

**Proof.** We consider the quadratic Lyapunov function

$$V(\mathbf{x}, \boldsymbol{\jmath}) = \sum\_{i=1}^{N} (\alpha\_i^2 + \boldsymbol{\jmath}\_i^2 / \alpha^2).$$

The derivative of V along the trajectories of the system (1) gives

$$\dot{V} = \sum\_{i=1}^{N} \left( 2\chi\_i \dot{\varkappa}\_i + 2\jmath\_i \dot{\wp}\_i / \alpha^2 \right) \tag{3}$$

which after substitution of the derivatives x˙<sup>i</sup> and y˙<sup>i</sup> , can be bounded from above as follows

$$\begin{split} &=\sum\_{i=1}^{N}\left(2x\_{i}y\_{i}-2y\_{i}\frac{\alpha^{2}x\_{i}}{\alpha^{2}}+2\frac{y\_{i}}{\alpha^{2}}\left(\lambda-x\_{i}^{2}-\frac{y\_{i}^{2}}{\alpha^{2}}\right)y\_{i}\right) \\ &\quad -2\mu x\_{i}\sum\_{j=1}^{N}a\_{ij}(x\_{i}-x\_{j})\right) \\ &\leq 2\frac{\lambda}{\alpha^{2}}\sum\_{i=1}^{N}y\_{i}^{2}-2\mu\sum\_{i,j=1}^{N}x\_{i}a\_{ij}(x\_{i}-x\_{j}) \\ &=2\frac{\lambda}{\alpha^{2}}\sum\_{i=1}^{N}y\_{i}^{2}-2\mu\sum\_{i,j=1}^{N}a\_{ij}(x\_{i}^{2}-x\_{j}x\_{i}) \\ &=2\frac{\lambda}{\alpha^{2}}\sum\_{i=1}^{N}y\_{i}^{2}-2\mu\sum\_{i,j=1}^{N}a\_{ij}(x\_{i}^{2}/2-x\_{j}x\_{i}+x\_{i}^{2}/2) \\ \end{split}$$

To show that the last term is negative definite, we use the balanced interconnection hypothesis. Let d<sup>j</sup> = P<sup>N</sup> i=1 aij = P<sup>N</sup> i=1 aji. Then,

$$\sum\_{i,j=1}^N a\_{ij} \mathbf{x}\_j = \sum\_{j=1}^N d\_j \mathbf{x}\_j = \sum\_{i=1}^N d\_i \mathbf{x}\_i = \sum\_{i,j=1}^N a\_{ij} \mathbf{x}\_i.$$

As a consequence,

$$\dot{V} \le \frac{\lambda}{\alpha^2} \sum\_{i=1}^N \wp\_i^2 - \mu \sum\_{i,j=1}^N a\_{ij} (\varkappa\_i - \varkappa\_j)^2 \le \frac{\lambda}{\alpha^2} \sum\_{i=1}^N \wp\_i^2.$$

Then the global part of statement follows by LaSalle invariance principle [40]. The local part follows by observing that for λ < 0 the linearization of model (1) at the origin is non-singular and therefore asymptotic stability of the origin implies that all eigenvalues have negative real part.

**Remark.** Because exponentially stability implies robustness to small perturbations, Theorem 2.1 remains true for small heterogeneity in the natural frequencies.

Note that, for λ < 0, the system in Equation (1) exhibits exponentially damped oscillations toward the origin (**Figure 1**, top panels).

Next we show that, at λ = 0 and identical natural frequencies model (1) undergoes a supercritical Hopf bifurcation inside the consensus space

$$\mathcal{C} = \{ (\mathbf{x}, \mathbf{y}) \in \mathbb{R}^{2N} : \mathbf{x}\_i = \mathbf{x}\_j, \ y\_i = y\_j, \ \forall i, j = 1, \dots, N \},$$

provided the graph is strongly connected. The linearization of the system (1) is given by

$$J = \begin{bmatrix} -\mu L & I\_N \\ -\alpha^2 I\_N & \lambda I\_N \end{bmatrix},\tag{4}$$

where I<sup>N</sup> is the N-dimensional identity matrix and L is the network Laplacian defined in section 2.1. Let 1<sup>N</sup> be the Ndimensional vector of ones. Given a (complex) vector ν = (v,w) in the consensus space C, that is, v = a1<sup>N</sup> and w = b1<sup>N</sup> for some a, b ∈ C, the eigenvalue problem for the Jacobian matrix Equation (4), restricted to the consensus space, takes the form

$$\begin{aligned} f\nu &= f \begin{bmatrix} a\mathbf{1}\_N \\ b\mathbf{1}\_N \end{bmatrix} \\ &= \begin{bmatrix} -\mu La\mathbf{1}\_N + b\mathbf{1}\_N \\ -\omega^2 a + \lambda b\mathbf{1}\_N \end{bmatrix} \\ &= \begin{bmatrix} b\mathbf{1}\_N \\ (-\omega^2 a + \lambda b)\mathbf{1}\_N \end{bmatrix} \\ &= \xi \begin{bmatrix} a\mathbf{1}\_N \\ b\mathbf{1}\_N \end{bmatrix} \end{aligned}$$

where ξ is a (complex) eigenvalue. In the third equality we used the fact that 1<sup>N</sup> is a right null eigenvector of L. The last equality shows that ν ∈ C is an eigenvector of J with eigenvalue ξ if and only if its components a and b satisfy

$$
\begin{bmatrix} 0 & 1 \\ -\alpha^2 & \lambda \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \xi \begin{bmatrix} a \\ b \end{bmatrix}. \tag{5}
$$

We can now easily solve Equation (5) to obtain the eigenvalues/eigenvectors pairs

$$\xi\_{\pm}(\lambda) = \frac{\lambda}{2} \pm \frac{1}{2}\sqrt{\lambda^2 - 4\alpha^2}, \quad \begin{bmatrix} a\_{\pm} \\ b\_{\pm} \end{bmatrix} = \begin{bmatrix} \frac{1}{\alpha^2}(\lambda - \xi\_{\pm}) \\ 1 \end{bmatrix} \tag{6}$$

For λ = 0 we got two purely imaginary eigenvalues which correspond to a supercritical Hopf bifurcation of model (1) inside the consensus space, as summarized in the following theorem and illustrated in (**Figure 1**, bottom panels).

**Theorem 2.2.** For almost all balanced, strongly connected interconnection topologies the following holds. For all µ > 0, the system (1) undergoes a supercritical Hopf bifurcation at λ = 0 with center manifold given by the consensus space C. Moreover, the family of periodic solutions born at the Hopf bifurcation are exponentially asymptotically stable and correspond to synchronous oscillations of the oscillator network.

**Proof.** By Theorem 2.1 the origin is locally exponentially stable for λ < 0. We further observe that, if the interconnection topology is strongly connected, then zero is a simple eigenvalue of L and therefore no either eigenvalue of J satisfies the same eigenvalue problem defined by Equation (5). It then follows by the center manifold theorem [41] and Equation (6), that the system (1) possesses a two-dimensional center manifold W<sup>c</sup> that is tangent to the consensus space C, for λ = 0. Moreover, this center manifold is exponentially attractive. By the Hopf bifurcation theorem [42], Equation (6) also implies that the system Equation (1) undergoes a supercritical Hopf bifurcation inside W<sup>c</sup> when λ crosses zero from negative to positive. By direct substitution inside the model equations, we see that along a generic member of the family of periodic orbits born at the Hopf bifurcation, oscillators are synchronously oscillating with each oscillator orbit given by (2).

Remark 2.3. Because Hopf bifurcation is codimension-zero (in the sense of [43]), it is persistent under small perturbations, which ensures that Theorem 2.2 remains true for small heterogeneity in the natural frequencies.

#### 2.2. Oscillator Death and Multi-Stability for Stronger Coupling

We now explore the phenomenon of "oscillator death," induced by strong coupling in model (1). We restrict our attention to the all-to-all coupling case, i.e., aij = 1 for all i 6= j. For λ > 0 and µ sufficiently small the synchronous oscillations born at Hopf bifurcation (Theorem 2.2) attract all trajectories. However, the system is multistable, as can be noted from the fact that increasing µ leads to the appearance of a family of steady states that attract some of the trajectories, but the synchronous periodic orbits remain locally exponentially stable (**Figure 2**). Indeed, depending on the initial conditions, only some trajectories converge to the synchronous oscillations. In the following we will provide geometric insights, without formal proof, about the mechanisms underlying oscillator death and multi-stability in model (1).

We start by observing that the oscillator death state is characterized by the presence of two dead oscillator clusters. Inside each cluster, oscillators converge to the same steady state. To analyze the appearance of oscillator death steady-states, we can simplify the model by assuming that (x<sup>i</sup> , yi) = (x1, y1) for all i = 1, . . . , N<sup>1</sup> and (x<sup>i</sup> , yi) = (x2, y2) for all i = 1, . . . , N2, where N1, N<sup>2</sup> < N, N<sup>1</sup> +N<sup>2</sup> = N, are the cluster sizes. The pairs (x1, y1) and (x2, y2) define the cluster states.

The cluster state dynamics can be easily derived and read

$$
\dot{\mathbf{x}}\_1 = \mathbf{y}\_1 + \mu \mathbf{N}\_2 (\mathbf{x}\_2 - \mathbf{x}\_1),
\tag{7a}
$$

$$\dot{\boldsymbol{\nu}}\_{1} = -\boldsymbol{\omega}^{2}\mathbf{x}\_{1} + \left(\boldsymbol{\lambda} - \boldsymbol{\kappa}\_{1}^{2} - \frac{\boldsymbol{\nu}\_{1}^{2}}{\boldsymbol{\alpha}^{2}}\right)\boldsymbol{\nu}\_{1},\tag{7b}$$

$$
\dot{\mathbf{x}}\_2 = \mathbf{y}\_2 + \mu N\_1 (\mathbf{x}\_1 - \mathbf{x}\_2),
\tag{7c}
$$

$$
\dot{\nu}\_2 = -\alpha^2 \varkappa\_2 + \left(\lambda - \varkappa\_2^2 - \frac{\wp\_2^2}{\alpha^2}\right) \wp\_2. \tag{7d}
$$

Each cluster state dynamics has the form

$$
\dot{\mathbf{x}} = \mathbf{y} + \mu \mathbf{N}\_{\dot{\mathbf{j}}} (\mathbf{x}\_{\dot{\mathbf{j}}} - \mathbf{x}),
\tag{8a}
$$

$$
\dot{\jmath} = -\omega^2 \varkappa + \left(\lambda - \varkappa^2 - \frac{\wp^2}{\alpha^2}\right) \wp,\tag{8b}
$$

where N<sup>j</sup> is the other cluster size and x<sup>j</sup> the other cluster state. A sufficient condition for the appearance of multiple steady states is that there must exist values of x<sup>j</sup> for which the model(8) has multiple steady-states. This condition can easily be verified by analyzing the dependence of the intersection of the nullclines of model (8) as a function of x<sup>j</sup> and the parameters µ and λ (**Figure 3**).

If the coupling strength is too small the origin is the only steady state (**Figure 3**). This steady state is unstable and all trajectories are attracted toward the synchronous periodic orbit. However, new steady states appear for larger values of µ. The critical value of µ for which the new steady-states appear can be found by computing the slope of the nullclines at the origin. The slope of the x-nullcline is evidently µN<sup>j</sup> . The slope of the xnullcline can be computed by implicit differentiation and is given by <sup>ω</sup> 2 λ . Multiple steady-states appear if µ > <sup>ω</sup> 2 Njλ (**Figure 2**).

#### 3. SYNCHRONIZATION AND FREQUENCY MODULATION IN LINEARLY COUPLED OSCILLATORS

In this section we present an alternative approach to study synchronization under the influence of noise using the Fokker-Planck equation (FPE). The modeling in this section can be thought of in the context of interacting populations of oscillators. Related work in the context of populations of synchronized neurons can be found in the work by Jiao et al. [44]. Introducing noise in the equation is natural from the biological perspective. However, even in the absence of noise, the introduction of random perturbations allows the extraction of information about the deterministic system. This is done by letting the perturbation amplitudes go to zero. To investigate the dependence of the synchronization frequency on the frequencies of the coupled oscillators, and the different coupling parameters, we assume

synchronization, which reduces the FPE equation to an equation in two variables.

We divide this section in two parts. First, we consider a simple deterministic system in which the effect of coupling can be understood. In the second part, we randomly perturb a more general version of the previous model to show that the FP equation provides an approximation for the syncrhonization frequency, and obtain some insights on the effect of noise. Let us then consider a system similar to the one already studied

$$\ddot{\boldsymbol{x}}\_{i} = -\boldsymbol{\omega}^{2}\boldsymbol{\chi}\_{i} + \nu \left[1 - \left(\boldsymbol{\chi}\_{i}^{2} + \frac{\dot{\boldsymbol{x}}\_{i}^{2}}{\boldsymbol{\omega}^{2}}\right)\right]\dot{\boldsymbol{x}}\_{i}$$

$$+\mu \sum\_{j=1}^{N} a\_{ij}\boldsymbol{\chi}\_{j}, \quad i = 1, \ldots, N,\tag{9}$$

Notice that x(t) = sin(ωt) is still a solution for the uncoupled system (µ = 0), independently of the value of ν. This system has the advantage of allowing direct calculations around the limit cycle, which can be written explicitly.

If we linearize the equation and take ν << 1, we can neglect the contribution of the disipative term. The resulting linear system is

$$
\ddot{\boldsymbol{\chi}}\_{i} = - (\boldsymbol{\omega}^{2} - \mu \boldsymbol{A}\_{i}) \boldsymbol{\chi}\_{i},
$$

where A<sup>i</sup> = P<sup>N</sup> j=1 aijx<sup>j</sup> , for i = 1, . . . , N. If all the aij = 1, corresponding to a fully connected network of oscillators, the previous system can be reduced to the equation

$$\ddot{\boldsymbol{x}} = -\left(\boldsymbol{\phi}^2 - \mu(\boldsymbol{N} - 1)\right)\boldsymbol{x},$$

assuming that a synchronized regime is established. This assumption may not always be biologically realistic, but it allows us to obtain the common synchronization frequency in a simple way. Later on we consider the general case and recover this formula as a particular case. This provides an estimate for the synchronization frequency of

$$
\Omega\_{sync} = \sqrt{\sigma^2 - \mu(N-1)}.
$$

Moreover, this reduction suggests that synchronized oscillatory behavior takes place for sufficiently small µ. That is, when

$$
\alpha^2 - \mu(N - 1) > 0\_\* 
$$

Otherwise, exponentially large growth can be expected. Notice that unless the aij are equal, the previous reasoning is not consistent and no conclusion can be drawn. We claim that introducing random perturbations and using the FP equation allows us to circumvent this difficulty and analyze the general case. This is the content of what follows. First of all, we write the system in the form

$$\begin{aligned} \dot{x}\_i &= y\_i \\ \dot{y}\_i &= f\_i(x\_i, y\_i) + \mu \sum\_{j=1}^N a\_{ij} x\_j, \quad i = 1, \dots, N. \end{aligned}$$

Notice that in the linearized regime, anologous to the reasoning for small ν in the previous example, we might naturally assume that

$$f\_i(x\_i, y\_i) \approx -a\_i^2.$$

Perturbing the equation with Brownian noise we have

$$\begin{aligned} dx\_i &= y\_i \, dt + \sqrt{2\varepsilon} \, dW\_{i1} \\ dy\_i &= \left[ -\omega\_i^2 x\_i + \mu \sum\_{j=1}^N a\_{ij} x\_j \right] dt + \sqrt{2\varepsilon} \, dW\_{i2}, \quad i = 1, \dots, N, \end{aligned}$$

where the Wij are uncorrelated Brownian motions for i, j ∈ {1, ..., N}. The probability density, u(x1, ..., xn, y1, ..., yn, t) of the system being in the state x1, ..., xn, y1, ..., y<sup>n</sup> at time t satisfies the FP Equation

$$\frac{\partial \mu}{\partial t} = \varepsilon \Delta \mu + \nabla (F(\varkappa, \jmath)\mu), \tag{10}$$

where F is the vector field determined by the right hand side of the stochastic system. Looking for stationary solutions, i.e., u<sup>t</sup> = 0 and explicitly substituting F in terms of x and y, the equation becomes

$$\left(\varepsilon \Delta \mu + \sum\_{i} \left(\wp\_{i} u\_{\mathbf{x}\_{i}} + (-\omega\_{i}^{2} \mathbf{x}\_{i} + \mu \sum\_{j \neq i} a\_{ij} \mathbf{x}\_{j}) u\_{\mathbf{y}\_{i}}\right)\right) = 0. \tag{11}$$

If we use the synchronization condition x<sup>1</sup> = ... = x<sup>n</sup> and y<sup>1</sup> = ... = yn, we obtain the equation

$$
\varepsilon \Delta + \eta \mu\_{\mathfrak{x}} + (-\sum\_{i} \alpha\_{i}^{2} + \mu \sum\_{i,j} a\_{ij}) \chi \mu\_{\mathfrak{y}} = 0.
$$

If we let

$$b = \sum\_{i} a\_i^2 - \mu \sum\_{i,j} a\_{ij,j}$$

we can write the Equation (11) as

$$
\varepsilon \Delta + n y u\_{\mathfrak{x}} - b \mathfrak{x} u\_{\mathfrak{y}} = 0.
$$

Assuming ε is small, it is reasonable to expect that the probability u will concentrate around the characteristic curves of the first order equation

$$a\gamma u\_{\mathfrak{x}} - b\mathfrak{x}u\_{\mathfrak{y}} = 0,$$

that are solutions to the system

$$\begin{aligned} \dot{\mathbf{x}} &= \eta \mathbf{y}, \\ \dot{\eta}\_i &= b \mathbf{x}. \end{aligned}$$

By taking the scalar product with the vector (x/n, y/b) we obtain the relation

$$
\dot{x}\frac{x}{n} + \dot{y}\frac{y}{b} = 0,
$$

or equivalently,

$$\frac{x^2}{n} + \frac{y^2}{b} = constant,$$

which defines the characteristic curves as ellipses. In turn, interpreting these as curves in the phase portrait, the resulting solutions would correspond to periodic trajectories with frequency

$$
\Omega\_{sync}^2 = -\frac{b}{n} = \frac{(\sum\_i \omega\_i^2 - \mu \sum\_{i,j} a\_{ij})}{n},\tag{12}
$$

which provides an estimate for the synchronization frequency in terms of the original frequencies and the coupling

parameters. Importantly, it shows that the synchronization frequency decreases with the coupling strength. In particular, formula Equation (12) can be used to study how coupled ultradian oscillations can give rise to circadian oscillations (**Figure 4**).

## 4. DISCUSSION AND SUMMARY

We have described, through basic geometrical analysis, the relationship between the dissipation coefficient, an intrinsic property of the oscillators we study, and the coupling strength µ in a network of diffusively coupled non-linear oscillators. Our analysis predicts the emergence of sustained oscillations for increasing values of the parameter λ in the system, only for a limited range of coupling strengths. It is reasonable to conjecture from this result, that there is a functional limit in the coupling strength for oscillating tissues in nature above which the tissue oscillations dies. To the best of our knowledge, it is the first time that diffusive coupling has been shown to be able to induce such oscillator death.

We have also derived an estimation for the synchronization frequency of a linearly coupled network of non-linear oscillators in terms of the oscillator natural frequencies and the coupling parameters [Equation (12)]. The presented results are indeed local, that is, the synchronous oscillation is only locally asymptotically stable. The oscillators are not synchronized at the beginning of the simulations, but their spread in state space is very small to ensure convergence to the synchronous oscillation. For a larger spread, a more complex behavior is observed.

We believe that these results constitute predictions that, although possibly difficult to test experimentally, would be worth verifying in light of the existing evidence about the joint frequency modulation of activity between different tissues during the day [23, 45].

The results we have presented thus far emphasize the importance of simple mathematical models in understanding situations where synchronization of multiple oscillating populations appears. The results presented here may help to shed light on both physiological and pathological phenomena involving synchronization of oscillators in different tissues (Parkinson's disease [46, 47], epilepsy [48, 49]). The other way around, it is also of potential importance to unravel mechanisms underlying the disappearance of coordinated oscillatory regimes. In a future publication, we plan to formally justify our estimations, and further, integrate the analysis of oscillations in the cellular and network levels of biological organization, to build up on our understanding of coupling oscillators at the tissue level. Two important extensions of the current models that we are studying are, the full characterization in higher codimension of the bifurcation structures of the system (1), and also, replacements of the van der Pol dynamics with biophysical models of excitable cells [50]. This last extension may prove useful to explain possible compensatory mechanisms that take place during the beginning of a pathology [51].

## AUTHOR CONTRIBUTIONS

All authors wrote the paper. ML-A proposed the original set of equations, AF and PP-L proposed revised sets of equations and extended models. All authors performed the analysis. AF and MH-V performed simulations.

## FUNDING

This research was supported by DGAPA-PAPIIT (UNAM) grants IA105518 and IA208618

## ACKNOWLEDGMENTS

The authors would like to thank Beatriz Fuentes-Pardo for her input in discussions about the modeling and its connections with the experimental work she has done. We would also like to thank Carolina Barriga-Montoya for her work on the Fokker-Planck equations and her support with the modeling.

#### REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Franci, Herrera-Valdez, Lara-Aparicio and Padilla-Longoria. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Stochastic Phylogenetic Algorithm for Mitochondrial DNA Analysis

M. Corona-Ruiz <sup>1</sup> , Francisco Hernandez-Cabrera<sup>1</sup> \*, José Roberto Cantú-González <sup>2</sup> , O. González-Amezcua<sup>1</sup> and Francisco Javier Almaguer <sup>1</sup> \*

<sup>1</sup> Facultad de Ciencias Físico-Matemáticas, Universidad Autónoma de Nuevo León, San Nicolás de los Garza, Mexico, <sup>2</sup> Escuela de Sistemas PMRV, Unidad Acuña, Universidad Autónoma de Coahuila, Saltillo, Mexico

This paper presents an exploratory analysis of the mitochondrial DNA (mtDNA) of 32 species in the subphylum Vertebrata, divided in 7 taxonomic classes. Multiple stochastic parameters, such as the Hurst and detrended fluctuation analysis (DFA) exponents, Shannon entropy, and Chargaff ratio are computed for each DNA sequence. The biological interpretation of these parameters leads to defining a triplet of novel indices. These new functions incorporate the long-range correlations, the probability of occurrence of nucleic bases, and the ratio of pyrimidines-to-purines. Results suggest that relevant regions in mtDNA can be located using the proposed indices. Furthermore, early results from clustering algorithms indicate that the indices introduced might be useful in phylogenetic studies.

#### Edited by:

Olcay Akman, Illinois State University, United States

#### Reviewed by:

Kyle B. Gustafson, Naval Sea Systems Command (NAVSEA), United States Monika Heiner, Brandenburg University of Technology Cottbus-Senftenberg, Germany

#### \*Correspondence:

Francisco Hernandez-Cabrera francisco.hernandezcbr@uanl.edu.mx Francisco Javier Almaguer francisco.almaguermrt@uanl.edu.mx

#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Genetics

Received: 01 April 2018 Accepted: 28 January 2019 Published: 08 March 2019

#### Citation:

Corona-Ruiz M, Hernandez-Cabrera F, Cantú-González JR, González-Amezcua O and Javier Almaguer F (2019) A Stochastic Phylogenetic Algorithm for Mitochondrial DNA Analysis. Front. Genet. 10:66. doi: 10.3389/fgene.2019.00066 Keywords: DNA, random-walk, Hurst exponent, detrended fluctuation analysis, Shannon entropy, coefficient of disequilibrium

## 1. INTRODUCTION

Previous mathematical studies on DNA sequences have seen a variety of approaches and frequently involve a numerical representation of the nucleotide chains. For instance, distance matrices have been constructed using different metrics (Randi et al., 2003; Liao and Wang, 2004; Zhang and Tan, 2007; Kandiah and Shepelyansky, 2013). These matrices, in combination with clustering methods, are used to evaluate phylogenetic relationships among species (Yu and Huang, 2013).

Other studies involve the representation of DNA sequences as random-walks, known as DNAwalks (Peng et al., 1994). The main objectives of these studies focus on the long-range correlations among nucleotides; i.e., "how the frequency of each nucleotide of a pairing nucleotide couple changes locally" (Namazi and Kiminezhadmalaie, 2015). These DNA-walk studies find differences in the long-range correlation between coding and non-coding DNA sequences (Peng et al., 1994).

Recently, DNA-walk analysis has been used in combination with the fractal dimension and Hurst exponent to identify mosaic structures in DNA that allow distinguishing between healthy and cancerous cells (Namazi and Kiminezhadmalaie, 2015).

Additionally, alternative statistical tools frequently used in DNA sequence analysis include Shannon entropy, which is a measure of the amount of "information" stored within a system (López-Ruiz et al., 1995). In a biological sense, Shannon entropy evaluates the probability of independent occurrences of each nucleic base in a DNA sequence. In recent studies, fluctuations in local Shannon entropy in DNA sequences have been analyzed to identify regions of repeating patterns of one or more nucleotides, known as tandem repeats (Thanos et al., 2018). The capability of Shannon entropy to highlight important segments in DNA sequences has led to the supported notion that entropy studies might be used for biological classifications of species (Melnik and Usatenko, 2014).

Similarly, the concept of complexity has played a central role in various DNA sequence analyses. For instance, López-Mancini-Calbet (LMC) complexity, employed in this paper, has led to the development of an effective gene-predicting technique (López-Ruiz et al., 1995; Monge and Crespo, 2015). In a recent study, the symbolic complexity of DNA sequences is used to identify segments resulting from random duplication, as well as changes in the speed of accumulation of point mutations (Salgado-Garcia and Ugalde, 2016).

Our objective is to examine the parameters previously mentioned to determine a small number of coefficients with biological relevance that may be used to determine rates of change in nucleotide bases, establish comparisons between regions, and better understand the relation among species in a phylogenetic sense.

This paper is structured as follows: section 2 introduces the concepts and methodology; section 3 presents the results obtained and the variables introduced; and section 4 is devoted to a discussion of the results, comments on the methodology in general, and final remarks. Tables and figures are incorporated in sections 2 and 3, respectively. The **Supplementary Material** includes a table with the identification codes for the data.

#### 2. METHODOLOGY

GenBank <sup>R</sup> is the National Institutes of Health's genetic sequence database made possible by the collaboration of several organizations. All datasets used within this work were obtained through GenBank because of its availability of access, encouragement of use, and the advantage that the information stays up-to-date.

A total of 32 complete mtDNA sequences of different species in the subphylum Vertebrata were selected. The lengths vary from 16, 207 to 18, 254 base pairs (bp). The choice of this type of DNA presents multiple advantages: it is relatively small in size (in contrast, human chromosomal DNA contains hundreds of millions bp); the sequences contain conserved regions, can be compared in blocks among different species, and contain a small percentage of non-coding regions; and the interpretation of the mutations in mtDNA as an estimator of evolutionary change (Barton and Jones, 1983). For these reasons, the exploratory nature of this study does not require additional information on the species themselves. Thus, the selection criteria focused on 32 different members from 7 groups intuitively related in taxonomic classes. The 32 NCBI codes from the data files have been attached in the **Table S1**.

A pre-processing of the data files consists of a realignment of the sequences to set the control region of the heavy chain (Hchain) in the direction of transcription as the new ending point. This realignment is done once. The displacement loop, or Dloop, is within the control region and the most varying region in mtDNA, with substantial differences observed even among individuals of the same species (Yamamoto, 2001). See **Figure S1** (Supplementary Material). Additionally, the header information was removed, which contains the identification key and the name of the organism. The downloaded files (in .fasta format) were processed using the programming language R version 3.4.4 (2018-03-15). The packages used are stringr and fractal.

## 2.1. DNA-Walk

DNA consists of sequences of nitrogenous bases: adenine (A), guanine (G), thymine (T), and cytosine (C). The length and distribution of the bases fluctuate from species to species. Several mappings have been introduced based on properties intrinsic to DNA. Moreover, adenine and guanine have a two-ring structure and belong to the purine group, while cytosine and thymine have a one-ring structure and belong to the pyrimidine group. Furthermore, adenine bonds with thymine through a double hydrogen bond, which is called a weak bond, while guanine and cytosine bond through a triple hydrogen bond, which is called a strong bond. **Figure 1** illustrates these descriptions. In summary, we have:


Considering the properties described previously, it is possible to read a DNA sequence and assign either a +1 or −1 depending on whether the respective nucleotide is a purine or pyrimidine (RY rule). This can be interpreted as random steps x<sup>i</sup> of a onedimensional walk. Then, the final position after n steps is given by

$$X\_{\mathfrak{n}} = \mathfrak{x}\_0 + \sum\_{i=1}^n \mathfrak{x}\_i \tag{1}$$

where x<sup>0</sup> = 0 by definition.

Let S = {s1s<sup>2</sup> . . .sM} be a nucleotide sequence of length M, where s<sup>k</sup> ∈ {A, C,G, T} for k ∈ {1, 2, . . . , M}. Hence, a one-dimensional DNA-walk can be defined through the following rules:

• RY rule:

$$x\_k = \begin{cases} 1 & \text{if } s\_k \in R = \{A, G\} \\ -1 & \text{if } s\_k \in Y = \{C, T\} \end{cases} \tag{2}$$

• SW rule:

$$\alpha\_k = \begin{cases} 1 & \text{if } s\_k \in \mathcal{S} = \{C, G\} \\ -1 & \text{if } s\_k \in \mathcal{W} = \{A, T\} \end{cases} \tag{3}$$

• KM rule:

$$\mathfrak{x}\_{k} = \begin{cases} 1 & \text{if } s\_{k} \in M = \{A, C\} \\ -1 & \text{if } s\_{k} \in K = \{G, T\} \end{cases} \tag{4}$$

where s<sup>k</sup> is the k−th nucleotide and x<sup>k</sup> is the value of the k−th assigned step in a DNA sequence. The path of the DNA-walk after n steps is then defined as the partial sums X<sup>n</sup> = x<sup>0</sup> + P<sup>n</sup> k=1 xk , where n ∈ {1, 2, . . . , M} and x<sup>0</sup> = 0.

In the context of DNA-walks, Equation (2) evaluates the tendency of changes between purines and pyrimidines. Transversions (substitutions of purines for pyrimidines, or vice versa) are less likely to happen and have been used to evaluate

molecular evolution (Stoltzfus and Norris, 2016). Thus, using this rule within corresponding blocks of nucleotides in different species, it is possible to observe changes in the DNA-walk that could be interpreted as an evolutionary variation. Similarly, Equation (4) is associated with the rate of recombination between transversions and transitions (purine-purine or pyrimidinepyrimidine substitutions).

Moreover, Equation (3) refers to the difference in abundance of the GC bond with respect to the AT bond. A higher GC content suggests a significantly higher temperature for DNA denaturing (melting temperature Tm). Previous studies have shown that GC content is associated to an age-related natural selection and environmental factors (Min and Hickey, 2008). Finally, it is assumed that each DNA-walk is an ergodic stochastic process. Specifically, the conceived notion adopted is that each DNA sequence may be used to represent the ensemble of DNA sequences of individuals within the same species.

In summary, the three assignment rules provide insight into the evolutionary aspects of the organisms considered.

#### 2.2. Hurst Exponent and DFA Exponent

Additional information of the long-range correlations of DNA-walks can be obtained via stochastic methods such as rescaled-range analysis and detrended fluctuation analysis. With these methods, it is possible to obtain the Hurst exponent, which represents a quantitative measure of the fractal nature of DNA sequences.

The Hurst exponent, here denoted by α, satisfies 0 < α < 1. In comparisons of mtDNA sequences, each Hurst exponent can be interpreted as a measure of the tendency of changes between nucleotides according to the rules mentioned in the previous section. The calculations used to obtain the Hurst exponent have been reported in previous studies (Peng et al., 1994; Buldyrev et al., 1995).

The Hurst exponent is directly related to the fractal dimension α ′ by the relation:

$$
\alpha' = 2 - \alpha. \tag{5}
$$

The fractal dimension evaluates changes in detail of the pattern of a DNA-walk with respect to the scale used for measurement.

An alternative method to calculate the Hurst exponent of a DNA-walk is DFA. In contrast to the rescaled-range analysis, DFA analyzes the random fluctuations of the DNA-walk without trend in the data (Peng et al., 1994; Buldyrev et al., 1995). The DFA exponent is computed using the following algorithm:

• Given a numerical sequence X = {X1, X2, . . . , XM}, calculate the cumulative sum

$$\mathcal{y}\_k = \sum\_{i=1}^k (\mathbf{X}\_i - \overline{\mathbf{X}}) \tag{6}$$

where k = 1, 2, . . . , M and X is the mean value of X.


$$F^2(L) = \frac{1}{M} \sum\_{k=1}^{M} \left| \mathcal{y}\_k - \mathcal{y}\_{kL} \right|^2. \tag{7}$$

• The slope β of the linear regression analysis in the scale log F(L)/ log L is an estimator of the Hurst exponent.

This method tests for self-similarity at different window sizes L. No correlation (or short-range correlations) gives stochastic properties such as those of a random-walk, so β = 0.5; in contrast, long-range correlations give a value of β 6= 0.5. Specifically, correlation yields β > 0.5, while anti-correlation gives β < 0.5.

This paper adopts a minimum block size of 4 nucleotides, while the maximum is B = M 2 , corresponding to half the length of the sequence in question. Should M be odd, B is rounded down.

#### 2.3. Chargaff Ratio

In a remarkable discovery, Erwin Chargaff determined that there is a balance held in DNA by the nucleobases (Chargaff, 1950), known as Chargaff's Rule. These state: (1) that globally (i.e., considering both strands of DNA) adenine is equal to thymine in quantity, and (2) that guanine is equal to cytosine in quantity. This result was the basis for the Watson-Crick model, which determined that adenine binds with thymine and that guanine binds with cytosine (Watson and Crick, 1953).

On this basis, and in the context of this work, the Chargaff ratio is defined as the ratio of pyrimidines to purines:

$$\xi = \frac{N\_C + N\_T}{N\_A + N\_G} \tag{8}$$

where NC, NT, NA, N<sup>G</sup> represent the amount of cytosine, thymine, adenine, and guanine, respectively, within one strand of DNA. Note that this value is always positive. If 0 ≤ ξ < 1, there are more purines than pyrimidines (i.e., N<sup>C</sup> + N<sup>T</sup> < N<sup>A</sup> + NG); similarly, ξ > 1 reflects an excess of pyrimidines over purines. A Chargaff ratio with value 1 results from an equal number of either type of nucleotide bases.

#### 2.4. Shannon Entropy

In his seminal paper, Claude Shannon introduced the concept of information entropy. It measures the "amount" of information or uncertainty of a system (Shannon and Weaver, 1998). Let = {ω1, ω2, . . . , ωN} be a set of events where each ω<sup>i</sup> has probability of occurrence p<sup>i</sup> ∈ [0, 1], for i = 1, 2, . . . , N. Thus, the Shannon entropy of the system is defined as

$$\mathcal{H} = -K \sum\_{i=1}^{N} p\_i \log\_2(p\_i),\tag{9}$$

where K is a positive constant chosen appropriately according to the units desired for measurement (thus, for this work, K = 1). For the case when p<sup>i</sup> = 0, p<sup>i</sup> log<sup>2</sup> (pi) = 0 in the limit definition. Also, note that the logarithm is in base 2; this is because information in a computer is encoded in binary digits, or bits, which are the basic units of measurement of information.

For N = 2, events ω<sup>1</sup> and ω<sup>2</sup> have probability p and 1 − p, respectively, see **Figure S2** (Supplementary Material). Thus, it can be seen that a maximum is attained at p = 1 − p = 1 2 . This result can be extended to the general case with N events. The proof requires Jensen's inequality for a concave function (in this case, the logarithmic function), and is given below. Using some algebra to rewrite Equation (9) with K = 1 yields

$$\mathcal{H} = \log\_2\left(\prod\_{i=1}^N \left(\frac{1}{p\_i}\right)^{p\_i}\right)$$

By the weighted arithmetic-mean and geometric-mean inequality, this implies that

$$2^{\mathcal{H}} = \prod\_{i=1}^{N} \left(\frac{1}{p\_i}\right)^{p\_i} \le \sum\_{i=1}^{N} p\_i \left(\frac{1}{p\_i}\right) = N$$

where equality (the maximum) is satisfied when p<sup>1</sup> = p<sup>2</sup> = · · · = pN. That is, when

$$\mathcal{H} = \log\_2(N). \tag{10}$$

To evaluate Shannon entropy in the context of DNA sequence analysis, it seems rather reasonable to define the set of possible events as = {A,G, C, T}. However, it is expected that the probability of occurrence of each nucleotide in a DNA sequence will likely be different for different species; thus, these associated probabilities will be calculated empirically for each DNA sequence in a straightforward fashion. That is, by counting the amount of each nucleotide within the sequence and taking the corresponding proportion by dividing by the total amount of nucleotides M. Thus, the probabilities will be given by

$$\rho\_A = \frac{N\_A}{M}, \qquad \rho\_C = \frac{N\_C}{M}, \qquad \qquad \rho\_G = \frac{N\_G}{M}, \qquad \qquad \rho\_T = \frac{N\_T}{M}, \tag{11}$$

where NA, NC, NG, N<sup>T</sup> are the amount of adenine, cytosine, guanine, and thymine, respectively.

In the context of DNA sequence analysis, maximum entropy is attained whenever the nucleic bases within a DNA sequence are found with equiprobability. It may thus be interpreted that such a sequence is the result of a random combination of these events. Any departure from the maximum value of the Shannon entropy due to an underlying structure might contribute to determining any tendencies present in a sequence, see **Figure S3** (Supplementary Material).

In a more general sense, the entropy fluctuations could be analyzed by means of the Local Shannon entropy. By studying the local fluctuations of entropy at a given scale, and across scales, an "entropic microscope" could highlight areas with a high degree of variation or, equally interesting, low degree of variation, as seen in previous studies (Melnik and Usatenko, 2014; Thanos et al., 2018).

#### 2.5. Coefficient of Disequilibrium

Additional information of DNA sequences can be derived from the deviations from equiprobability of occurrence of each

FIGURE 2 | DNA-walk illustration for various species using the purine-pyrimidine rule. Observe the vicinity of nucleotide 2, 700 and the change in tendency from a purine-rich region (positive slope) to a predominance of pyrimidines for the remaining DNA-walk (negative slope).

nucleotide. This measure is known as disequilibrium (López-Ruiz et al., 1995). The events in the set have probability p<sup>i</sup> for i = 1, 2, 3, 4. The coefficient of disequilibrium, D, is defined as:

$$\mathcal{D} = \sum\_{i=1}^{N=4} \left( p\_i - \frac{1}{4} \right)^2. \tag{12}$$

FIGURE 3 | DNA-walk illustration for various species using the strong- and weak-bond rule. Observe the immediate (and consistent) tendency. This indicates that mtDNA is rich in adenine and thymine, whose type of bond is weaker than that of cytosine and guanine.

rule. The figure shows a higher amount of adenine and cytosine.

This sum of squared distances can be seen as a type of variance. Note that D = 0 in the case of equilibrium. Any deviation from this would result in D > 0. The maximum disequilibrium value, Dmax = 3 4 can be obtained using multivariate calculus.

The coefficient of disequilibrium may represent a measure of relatedness between a DNA sequence and one resulting from a random process if each (independent) event has a probability p<sup>i</sup> of occurrence. That is, larger deviations from an equiprobable space yield higher coefficients of disequilibrium. It can be observed that this behavior counters that of the Shannon entropy in an intuitive manner.

## 2.6. Coefficient of Complexity

The coefficient of complexity C is then given by the product of the Shannon entropy (9) and the coefficient of disequilibrium (12), as in (13). It can be seen from (12) that D resembles the definition of

TABLE 1 | Results of the Chargaff ratio and Shannon entropy for all groups.


TABLE 2 | Results of the Hurst exponent for all groups and each of the three random-walk rules.

TABLE 3 | Results of the DFA exponent for all groups and each of the three random-walk rules.



$$\mathcal{C} = \mathcal{H}\mathcal{D} = \left(-\sum\_{i=1}^{N} p\_i \log\_2(p\_i)\right) \left(\sum\_{i=1}^{N} \left(p\_i - \frac{1}{N}\right)^2\right). \tag{13}$$

variance; thus, the coefficient of complexity can be interpreted as a measure of dispersion within the information stored in a system (López-Ruiz et al., 1995).



#### 3. RESULTS

The three DNA-walks for the 7 groups are depicted in **Figures 2**– **4**. Results for the Chargaff ratio ξ and Shannon entropy H are shown in **Table 1**, while **Tables 2**, **3** contain the Hurst and DFA exponents for each type of random-walk and for each sequence.

In **Figure 2**, there is an initial upward trend that is present irrespective of the species. The RY rule (Equation 2) implies that a (local) inclination toward the positive direction of the vertical axis corresponds to a (local) majority of purines (adenine or guanine). Similarly, the downward trend in **Figure 3** reflects a consistent predominance of the weakly-pairing bases, adenine or thymine (considering rule SW). Thus, adenine dominates within the range 0− ∼ 3, 000 bp.

The Hurst exponents for the rules RY, SW, and KM (Equations 2–4, respectively) fall in the range of 0.900 − 0.912 and imply a long-term positive autocorrelation. To put it into perspective, a Hurst exponent value of 0.9 indicates that, on average, the tendency of changes between nucleotides varies slightly as the sub-sequence size is changed. Moreover, the proximity of the Hurst exponent toward unity suggests that either purines or pyrimidines are predominant; it cannot distinguish, however, which one prevails. Similarly, the DFA exponents fall within 0.64 − 0.91 which implies the existence of strong long-range correlations in the sequences even after detrending. Interestingly, neither the Hurst nor DFA exponent values are near zero in any of the species considered. A possible explanation is that the tendency of changes between nucleotides does not vary randomly; i.e., mtDNA has an informational structure.

For all the DNA sequences, the Chargaff ratio is positive with ξ > 1, implying a larger amount of pyrimidines than purines. This implication is visually reflected in the overall downward tendency of the curves in **Figure 2**.

The disequilibrium coefficient takes values D ∈ (0.01 − 0.03). From Equation (12), values near 0 imply that the probabilities p<sup>i</sup> ≈ 1 4 for any of the four nucleic bases. In other words, the disequilibrium values obtained suggest that the four nucleotide bases appear with almost the same proportion within each of the 32 mtDNA sequences. This is further supported by the Shannon entropy values. In this case, Equation (10) and N = 4 yield a (theoretical) maximum entropy value <sup>H</sup> = log<sup>2</sup> (4) = 2. Hence, the empirical entropy values H ∈ (1.89 − 1.97) suggest near-equiprobability among the nucleic bases.

A graph of D vs. the Shannon entropy H suggests a linear relation. On this account, the disequilibrium coefficient is omitted for the remainder of the study. In addition, the complexity coefficient is omitted due to its direct proportionality to D. See **Figure S4** (Supplementary Material).

This work proposes three new evolutionary indices as functions of Shannon entropy, the Chargaff ratio, and the fractal dimensions derived from the Hurst and DFA exponents:

$$\nu\_1 = \mathcal{H} \* \log \left[ \alpha'\_{RY} \* \xi \* \log \left( \alpha'\_{KM} \right) \right] \tag{14}$$

$$\nu\_2 = \log \left[ \boldsymbol{\beta}\_{RY}^{'} \* \log \left( \boldsymbol{\beta}\_{KM}^{'} \right) \right] \tag{15}$$

$$\nu\_3 = \log \left[ \boldsymbol{\beta}\_{SW}^{'} \* \log \left( \boldsymbol{\alpha}\_{SW}^{'} \right) \right]. \tag{16}$$

These indices reflect the long-range correlations found in DNAwalks and the information given by Shannon entropy and the Chargaff ratio.

The fractal dimensions α ′ and β ′ are derived from the Hurst and DFA exponents, respectively, using Equation (5).

The natural logarithm can be seen as a transformation that maximizes the differences between the coefficients. Equations (14), (15), and (16) are defined from an evolutionary perspective, while Equation (16) provides information on the energy content of sequences.

In Equation (14), the logarithm of the fractal dimension derived from the Hurst exponent using the KM rule provides information regarding the transversions and transitions of the entire DNA sequence. On the other hand, the Chargaff ratio is used as a weighting factor for the fractal dimension derived using the RY rule. The logarithm of the product of these quantities provides an evolutionary measure related to the long-range correlations. The last term in the equation (the Shannon entropy) evaluates the probability of independent nucleotide changes for a given DNA sequence.

Equation (15) uses the fractal dimensions of the DFA exponents, which are computed using the detrended DNAwalks. Therefore, it is not accurate to include the Chargaff ratio or Shannon entropy as normalization parameters. Finally, Equation (16) represents a measure of the natural selection factors in relation to the environment. Results for v1, v2, v<sup>3</sup> are shown in **Table 4**.

Clustering algorithms may benefit from the proposal. Preliminary results, shown in **Figure 5**, suggest a possible application in studies centering on the evolutionary relations among species. The proposed indices are used in the groupaverage agglomerative clustering algorithm with Euclidean metric and the sum of distances as the clustroid. Furthermore, an additional grouping was constructed using a traditional program, ClustalW, which is frequently applied to the study of phylogenetic trees, as seen in **Figure 6**.

The implementation of the algorithm using the R programming language is not computationally demanding, with running times of about 15–20 min. In comparison, ClustalW requires about 2 and a half hours for the construction of the phylogenetic tree of 32 mtDNA sequences.

The comparative analysis between the two methods shows consistency among the group of primates and other mammals sharing a common ancestry of similar lineage to the lemur. On the other hand, the marsupials and rodents (including the common rabbit) are more closely grouped with the stochastic algorithm and present a common ancestor, just as calculated by the traditional method. Other groups that share proximity with both methods are the reptiles and the birds, as well as the fish group and some amphibians.

The most pronounced differences are found in certain taxa. The proposed method relates the rabbit more closely to rodents, with characteristics similar to marsupials. Meanwhile, the traditional method positions the rabbit closer to primates. Another interesting point is that the proposed stochastic method shows that small reptiles and birds are more closely related, while the traditional method relates the birds closer to large reptiles.

## 4. CONCLUSIONS

As has been suggested by other studies, Shannon entropy and Hurst and DFA exponents provide insight into the properties of DNA sequences (Peng et al., 1994; Oiwa and Glazier, 2004; Melnik and Usatenko, 2014; Monge and Crespo, 2015; Namazi and Kiminezhadmalaie, 2015; Salgado-Garcia and Ugalde, 2016; Thanos et al., 2018). This exploratory analysis combines various

measures utilized in the literature to establish a biologically meaningful measure of distinction among species.

Our proposal defines new indices as functions of Shannon entropy, the Chargaff ratio, and fractal dimensions using rescaled-range analysis and DFA. These indices can be employed to construct phylogenetic trees using clustering algorithms.

Long-range correlations attributed to DNA-walks can be identified during our study. These can represent data with persistence in its evolutionary memory; i.e., that mtDNA sequences contain highly conserved regions among similar species.

The comparison between the traditional and the proposed clustering method shows clear agreements; however, there are differences that must be analyzed under an evolutionary perspective. For example, we notice that the mtDNA sequences of the common rabbit and the common snapping turtle show different properties in both methods. According to the established phylogeny, the placement of the rabbit is closer to the rodents. Interestingly, results of the stochastic hierarchical clustering suggest a potential application for phylogenetic studies.

Evolutionary processes are associated to an adaptive selection of the species throughout millions of years. However, the fluctuations of the changes in nucleotide bases could be random in order to find new sequence combinations. The proposed method attempts to measure the stochastic fluctuations to yield indices that allow the observation of tendencies and correlations in the mutations that produce new species throughout evolutionary history.

## REFERENCES

Barton, N., and Jones, J. (1983). Mitochondrial DNA: new clues about evolution. Nature 306, 317–318. doi: 10.1038/306317a0

## AUTHOR CONTRIBUTIONS

MC-R provided data collection of the mtDNA sequences from the GenBank <sup>R</sup> , worked on numerical and graphical results, and drafted the article. FJ provided numerical analysis, methodology, and mathematical insight. FH-C rendered numerical analysis, as well as mathematical and biological interpretations. JC-G contributed with numerical analysis, revision, critical revision for important intellectual content, and co-final approval of the version to be published. OG-A provided textual and structural revision of the co-final version of this work.

## ACKNOWLEDGMENTS

Thanks are due to the Consejo Nacional de Ciencia y Tecnología (Conacyt) for providing a scholarship for one of the authors. Special recognition is given to the Universidad Autónoma de Nuevo León, the Facultad de Ciencias Físico-Matemáticas, and the Centro de Investigación en Ciencias Físico-Matemáticas for logistical support given during our research endeavors. Thanks are due to Programa para el Desarrollo Profesional Docente, para el Tipo Superior (PRODEP) for the support for the publication of the article.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2019.00066/full#supplementary-material

Buldyrev, S. V., Goldberger, A. L., Havlin, S., Mantegna, R. N., Matsa, M. E., Peng, C.-K., et al. (1995). Long-range correlation properties of coding and noncoding dna sequences: genbank analysis. Phys. Rev. E 51, 5084–5091. doi: 10.1103/PhysRevE.51.5084


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Corona-Ruiz, Hernandez-Cabrera, Cantú-González, González-Amezcua and Javier Almaguer. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Probabilistic Assessment of Investment Options in Honey Value Chains in Lamu County, Kenya

Joshua Wafula<sup>1</sup> \*, Yusuf Karimjee<sup>1</sup> , Yvonne Tamba<sup>1</sup> , Geoffrey Malava<sup>1</sup> , Caroline Muchiri <sup>1</sup> , Grace Koech<sup>1</sup> , Jan De Leeuw1,2, Josephat Nyongesa<sup>1</sup> , Keith Shepherd<sup>1</sup> and Eike Luedeling1,3

<sup>1</sup> World Agroforestry Centre (ICRAF), Nairobi, Kenya, <sup>2</sup> ISRIC World Soil Information, Wageningen University, Wageningen, Netherlands, <sup>3</sup> Department of Horticultural Sciences, University of Bonn, Bonn, Germany

Designing and implementing biodiversity-based value chains can be a complex undertaking, especially in places where outcomes are uncertain and risks of project failure and cost overruns are high. We used the Stochastic Impact Evaluation (SIE) approach to guide the Intergovernmental Authority on Development (IGAD) on viable investment options in honey value chains, which the agency considered implementing as an economic incentive for communities along the Kenya-Somalia border to conserve biodiversity. The SIE approach allows for holistic analysis of project cost, benefit, and risk variables, including those with uncertain and missing information. It also identifies areas that pose critical uncertainties in the project. We started by conducting a baseline survey in Witu and Awer in Lamu County, Kenya. The aim of the survey was to establish the current farm income from beekeeping as a baseline, against which the prospective impacts of intervention options could be measured. We then developed an intervention decision model that was populated with all cost, benefit and risk variables relevant to beekeeping. After receiving training in making quantitative estimates, four subject-matter experts expressed their uncertainty about the proposed variables in the model by specifying probability distributions for them. We then used Monte Carlo simulation to project decision outcomes. We also identified variables that projected decision outcomes were most sensitive to, and we determined the value of information for each variable. The variable with the highest information value to the decision-maker in Witu was the honey price. In Awer, no additional information on any of the variables would change the recommendation to invest in honey value chains in the region. The analysis demonstrates a novel and comprehensive approach to decision-making for different stakeholders in a project where decision outcomes are uncertain.

Keywords: value chains, probabilistic projection, decision outcomes, uncertainity, value of information

## INTRODUCTION

How to approve and prioritize among projects that aim at biodiversity conservation has been highlighted as one of the most critical decisions that conservation planners face [1]. This is not surprising, because conservation outcomes are often achieved through complex mechanisms, and the success of conservation actions is rarely guaranteed, with many uncertainties preventing precise

#### Edited by:

José Roberto Cantú-González, Universidad Autónoma de Coahuila, Mexico

#### Reviewed by:

Maria Elena De Giuli, University of Pavia, Italy Barret Pengyuan Shao, Independent Researcher, United States

\*Correspondence: Joshua Wafula jokumu56@yahoo.com

#### Specialty section:

This article was submitted to Environmental Informatics, a section of the journal Frontiers in Applied Mathematics and Statistics

> Received: 15 December 2017 Accepted: 09 March 2018 Published: 27 March 2018

#### Citation:

Wafula J, Karimjee Y, Tamba Y, Malava G, Muchiri C, Koech G, De Leeuw J, Nyongesa J, Shepherd K and Luedeling E (2018) Probabilistic Assessment of Investment Options in Honey Value Chains in Lamu County, Kenya. Front. Appl. Math. Stat. 4:6. doi: 10.3389/fams.2018.00006 impact prediction. Success is even harder to predict, when conservation agencies aim to strengthen biodiversity indirectly, e.g., by supporting livelihoods and economies of local people as an incentive for them to conserve biodiversity outcomes [2]. Investing in biodiversity based value chains does not necessarily result in positive biodiversity outcomes. Negative impacts can arise, when value chain development results in depletion of the biodiversity that forms the resource base, on which the value chain depends (e.g. fisheries or non-timber forest products).

The production of honey is an example of a biodiversity based value chain that strengthens rather than erodes the conservation of biodiversity [3]. This is because honey producers have an interest to conserve the vegetation and plant species that produce the nectar and pollen that supports the value chain. The development of honey value chains typically revolves around a combination of introducing improved bee keeping and honey production techniques and improved access to markets for honey [4]. Yet, while attractive at first sight, such improved techniques are not always easily adopted [5]. An important reason for this is uncertainty among farmers about the financial outcomes of their investment in improved honey production techniques.

A detailed cost-benefit analysis on beekeeping projects can be considered to reduce the perceived uncertainty. However, there are rarely sufficient data on all relevant aspects of an investment decision to allow precise, purely data-driven projections to support decision-making [6]. Given such a lack of perfect knowledge, decision-makers need appropriate tools for handling uncertainties, and for identifying and prioritizing knowledge gaps, whose narrowing would reduce their chance of selecting a suboptimal decision option [7, 8]. Furthermore, decisionmakers need improved capabilities to quantify risks surrounding proposed interventions, because failure to adequately account for risk can lead to high chances of project failure [9].

The Stochastic Impact Evaluation (SIE) approach allows for a structured decision analysis that incorporates all relevant variables, even those with uncertain and missing information [10]. It considers risk factors that may compromise project success or affect project performance. The approach incorporates Value of Information analysis that prioritizes critical uncertainties in a project, where further research has the greatest potential of enhancing clarity on the decisions. The present study uses the SIE approach to assess investment decisions in honey value chains for the Intergovernmental Authority on Development (IGAD) in its program on Biodiversity Management (BMP).

#### Study Background

IGAD-BMP partnered with the World Agroforestry Centre (ICRAF) to implement the program's biodiversity-based interventions along the Kenya-Somalia border. During the project inception phase, stakeholders were consulted, and they proposed participatory honey value chain development as one of the economic incentives for biodiversity conservation. Communities from Witu and Awer were selected to pilot the beekeeping project in Lamu County. One hundred farmers from both communities were selected and trained in beekeeping techniques to boost their honey production knowledge and improve their access to formal markets.

The training was held in May 2015 in Malindi, Kenya. Malindi was selected for training due to its proximity to the Arabuko-Sokoke Forest, where surrounding communities intensively practice modern beekeeping. This provided an opportunity for the trainees to learn from established beekeepers. The training brought together different honey value chain actors, including representatives from the Ministry of Livestock and Fisheries in Kenya. During the workshop, stakeholders agreed that honey value chains had potential for being a viable investment option for Witu and Awer communities. A baseline survey (**Figure 1**) was commissioned to establish the current net income from beekeeping, the actual number of beekeepers among the trained farmers, beekeeping practices and risks associated with beekeeping.

In Awer, the average net income from beekeeping was 82 thousand Kenyan Shillings (Ksh) per year, approximately Ksh 225 per day (**Table 1**). Traditional beekeeping was widely practiced and characterized by high productivity and wellestablished markets. In Witu, the average net income from beekeeping per farmer was Ksh 39 thousand (approximately Ksh 107 per day). Very few farmers were currently practicing beekeeping. The region also had very low honey productivity and low farm gate prices, although a higher percentage of interviewed farmers were educated compared to Awer.

The baseline survey informed the project implementation approach for both communities. In Awer, it seemed sensible to support traditional beekeeping, since the region was characterized by high honey productivity and favorable farm gate prices. Establishment of a honey collection center and support of value addition activities, such as packaging, branding and processing of other bee products, were the most promising investment options. In Witu, the favored approach was to invest in modern beekeeping with emphasis on intensive farmer training and engagement of all honey value chain actors to boost productivity. Bridging existing gaps between beekeepers and formal markets was also a priority for the implementer to ensure favorable farm gate prices for the beekeepers.

While clear strategies for HVC development in the two regions thus emerged from the initial consultations and baseline survey, IGAD-BMP was still lacking certainty that the investment would raise farm incomes. To clarify the prospects for the intervention, the analysis aimed at (i) projecting impacts on beekeeping farms using a probabilistic impact evaluation approach, (ii) identifying and mitigating critical uncertainties in the project, and (iii) providing recommendations for project monitoring to reduce the chance of negative outcomes.

## MATERIALS AND METHODS

We used ICRAF's Stochastic Impact Evaluation approach (**Figure 2**), which is based on the principles of Applied Information Economics (AIE) [9, 11], to project the impact of the decision on different stakeholders in the project. The method makes use of all available information, including expert

TABLE 1 | Summary of the baseline survey results for net income from beekeeping in Witu and Awer.


knowledge, in making impact projections. The SIE approach can be applied even in the absence of perfect data. We quantified all risks, costs and benefits based on the current state of knowledge about them, considering the uncertainties about intervention outcomes that result from this. We also embraced the concept of Value of Information, which is useful for determining decisionspecific knowledge gaps that decision-supporting research should address [12].

## Modeling Process

The SIE process begins by thoroughly defining the decision to be made. According to Hubbard [9], a decision may merit structured decision analysis if it has two or more realistic alternatives, invokes some form of uncertainty or dilemma, has potentially negative consequences if it turns out that the wrong position was taken, and involves a decision-maker. For this study, HVC development was identified as a viable investment option during the consultative workshop, but stakeholders were still unsure that the decision to invest in HVCs would benefit them. This led to the question: "Should IGAD-BMP invest in the proposed honey value chain intervention as an economic incentive for Witu and Awer communities to conserve biodiversity?"

To clearly define this question, a decision analysis team was constituted with the aim of clarifying whether IGAD-BMP was interested in return on investment, who exactly the agency was targeting to raise incomes from beekeeping, and how the project was going to be implemented. The team that was convened to develop the HVC decision model consisted of 8 members: IGAD-BMP's HVC intervention project manager, the biodiversity management program coordinator, two decision analysis experts, one participatory modeling scientist and four research assistants. The initial conceptual model (**Table 2**), which aimed to reflect the overall structure of the decision at hand, was then sent to a principal bee health scientist from the International Centre of Insect Physiology and Ecology (icipe) for review. The updated model (**Figure 3**) was translated into a set of equations that represented the team's understanding of the decision's impact pathway [13].

Decision outcomes were projected using the decisionSupport package [14] for the R programming environment [15]. For this, the model equations were coded as a welfare function in R programming language (Data Sheets 3, 4). The model was parameterized using estimated probability distributions for variable values (Data Sheets 1, 2) (often specified as 90% confidence intervals), which were obtained from literature review and five established beekeepers from Mwingi in Kitui County, Kenya. The process required the team to undergo calibration training—a process that improves the capacity to make range estimates for which one is 90% confident that the actual value lies within the provided range [9]. This was done for all variables, including those for which no other information was available. Using the estimated variables as inputs, the decision model was run as a Monte Carlo simulation [16, 17]. This approach produced a distribution of possible decision outcomes by running the model a large number of times, each time fed with a different set of random draws from the defined distributions for the variables in the equations [18]). To identify variables that were most uncertain in the analysis, we used Partial Least Squares (PLS) regression [19], which was also implemented in the decisionSupport package [14].

The most uncertain variables in a project, which outcome projections are most sensitive to, are not necessarily most important to the decision-maker, because new information on them may not be able to change the recommendation emerging from the decision model. The value of additional information on a variable for decision-makers is determined by whether this information has the potential to change the sign of the expected value of the decision, which would change the preferred decision option [9].

To identify high-value variables, the Expected Value of Perfect Information (EVPI) was calculated by detecting correlations


between input and output variables and identifying those that were of importance to a decision-maker, i.e., the variables, whose measurement could help reduce the expected opportunity loss of the decision. Since there was no basis for deciding on the functional form of the relationship between test and outcome variable, a non-parametric test was used to detect these relationships. We used Spearman's rank correlation, using the usual statistical cutoff criterion of α < 0.05 to exclude all variables that were not correlated with projected decision outcomes. For those that were not correlated, the value of information was zero, since the variable had no significant relationship with the outcome. This criterion excluded most

variables, and discriminated well between informative and obviously uninformative variables.

We then sorted the dataset by values of each remaining input variable. In the sorted dataset, it was not possible to identify a clear threshold, where the expected decision outcomes transition from positive to negative. This was because the data was still too "noisy." To more clearly expose the effect of uncertain variables on the decision recommendation, we processed the dataset using a second-order low-pass Butterworth filter, with a critical frequency of one divided by one tenth of the number of values in the Monte Carlo output. This resulted in a smooth dataset, in which it was often possible to identify a threshold, where the sign of projected decision outcomes was reversed, e.g. for variable values above the threshold, positive outcomes would be expected, while values below the threshold would likely lead to negative outcomes (note that using a signal processing filter is a pragmatic way to solve the computational challenge of calculating the EVPI, but it introduces a small amount of inaccuracy into the threshold identification). With this threshold identified, the EVPI was computed as the sum of all outcomes with a sign that did not correspond to that of the expected value (e.g., all negative outcomes, when the analysis produced a positive expected value), multiplied by their respective chance of occurrence. This EVPI procedure was applied to all output variables.

## RESULTS

Due to the economic and biophysical disparities between Awer and Witu, the two regions were modeled separately. The outcome of the analysis was expressed as the net present value (NPV) of the intervention for a farmer and for the overall project. Emphasis was placed on the farmer NPV and overall project NPV, since the implementer's direct objective was to raise farmer incomes from beekeeping. Cash flows were also illustrated, and variables that had information value for the decision maker were identified.

## Witu Community

In Witu, the median of the modeled distribution of average annual monetary income per farmer practicing beekeeping was Ksh 65 thousand, with 90% confidence that the actual income lay within the range of Ksh −36 thousand to 140 thousand for a farmer who would continuously practice beekeeping (**Figure 4**). The chance of a negative average monetary NPV for a farmer in this region was 24%.

The honey price, the amount of honey produced per hive, the number of hives per farmer, the number of harvesting seasons per year and the honey processing cost were the most important variables according to the PLS analysis. The distribution of total farmer NPV, including indirect benefits from beekeeping, e.g. through improved pollination of crops, had a median of Ksh 72 thousand, with a 90% confidence interval of Ksh-30–251 thousand.

Since both positive and negative farmer outcomes appeared plausible according to the simulation results, the EVPI analysis indicated that additional information about four variables could potentially change the recommendation emerging from the SIE process. These were the honey price, with a value of information of about Ksh 2750, followed by honey production per hive (Ksh

1750), the number of hives per farmer (Ksh 1250) and the

indicate positive correlations of uncertain variables with the outcome variable, while red bars indicate negative correlations.

number of harvesting seasons per year (Ksh 50).

The distribution of the projected NPV for the project in Witu had a median of Ksh 36 million, with 90% confidence that the actual NPV for the project lay within the range of Ksh-14–150 million. The model responded most sensitively to the amount of honey produced per hive. Seven other variables also had important impact on projected outcome values (**Figure 5**).

The distribution of projected decision outcomes included both positive and negative outcomes, and additional information on some individual uncertain variables had potential to change the decision recommendation. EVPI analysis indicated that information on the honey price in Witu was the most valuable to the project implementers, with a value of about Ksh 1.1 million. The amount of honey produced per hive (Ksh 0.6 million) and the number of hives per farmer (Ksh 0.4 million) were also of value to the project implementers.

## Awer Community

In Awer, the distribution of projected average monetary income per farmer practicing beekeeping had a median of Ksh 130 thousand per year, with 90% confidence that the actual value lay within the range of Ksh 1.5–340 thousand per year. The chance of loss for a farmer in this region was 4.6%. Total farmer NPV including non-monetary benefits had a median of Ksh 140 thousand with a 90% confidence interval of Ksh 7.2–340 thousand. The price of honey was identified as the most uncertain variable by the PLS analysis, alongside five other variables (**Figure 6**).

EVPI analysis indicated that none of the variables had any information value for the farmer. This meant that, if the NPV was the main criterion for evaluating the attractiveness of the HVC intervention, then farmers would be well-advised to engage in beekeeping, and no additional information on any of the variables would change this recommendation.

The overall project NPV had a median of Ksh 38 million per year, with 90% confidence that the actual value lay within the range of Ksh-0.9–120 million. PLS analysis indicated that the honey price was the most uncertain variable in the projection of project NPV in Awer (**Figure 7**). Although PLS indicated that projected project outcomes responded to variation in a number of variables, no additional information on any of these variables had the potential to change the recommendation that the project implementer should invest in honey value chains in the region.

## DISCUSSION

Beekeeping is a high-risk, high-return venture that requires a well thought-out project design to maximize returns [20]. Introduction of the HVC intervention as an economic incentive for rural communities along the Kenya-Somalia border to conserve the environment may generate significant impact in terms of raising income (average of Ksh 65 thousand per farmer for Witu and Ksh 130 thousand per farmer for Awer). However, the success of this intervention, especially for regions such as Witu, requires close monitoring by the implementing agency to minimize the chance of loss (24%).

In Witu, the most critical variable was the price of honey, which individual beekeepers have relatively little influence on. Current honey price, based on the baseline survey, is within the range of 3–7 hundred Ksh per liter. At this price, returns may be too low for a farmer who pursues beekeeping as the only source of income. For a farmer who invests in beekeeping as a supplementary source of income, the current honey price looks promising. Analyzing the impact of different factors that may affect the honey price in the future, e.g., packaging, honey quality, access to formal markets and exploitation by middlemen, will provide useful information to support farmers' decisions on investments in honey value chains. The amount of honey produced per hive is also a critical variable to farmers in Witu (**Figure 4**). Since productivity is a function that depends on many

factors, such as good apiary management, access to extension services, technology deployed and bee forage availability, studies on these factors promise to reduce decision uncertainty for farmers. Modifications to the original project design that address the influence of these variables on the farmer NPV could reduce the chance of losses for beekeepers.

indicate positive correlations of uncertain variables with the outcome variable, while red bars indicate negative correlations.

Most model runs for project NPV in Witu indicate a positive outcome, although there is still a chance (23%) of not achieving the goal of raising incomes for farmers in the region. This is because beekeeping is a relatively new venture for most farmers. The value chain is not well developed and a lot of investment has to go into farmer training to increase productivity. The honey price, the amount of honey produced per hive and the number of hives per farmer are the most critical variables for the project implementer. It is therefore important for the implementer to work with farmers to acquire more information on these variables and how they will influence their respective NPVs.

In Awer, the HVC intervention can increase monetary income for beekeepers by about Ksh 48 thousand based on the baseline survey. This is slightly higher than in Witu (Ksh 26 thousand), but this is expected, since most farmers in the region were already practicing beekeeping. The project plan to support existing traditional beekeeping methods, establish a honey collection center and support value addition activities such as packaging, branding and processing of other bee products, was likely to be more effective in raising incomes, since it would focus on bridging gaps in the chain. In Witu, the focus is on introducing the HVC to the region.

The EVPI analysis indicated that no additional information on any of the variables could change the recommendation that farmers should invest in beekeeping (**Figure 6**). This is not surprising, since the region is characterized by high honey productivity and well-defined markets. The overall project NPV in Awer suggests that introducing HVC is economically

viable, even when fully considering the cost incurred by the implementing agency (**Figure 7**). The effort to improve certainty on the business case may not influence the decision to invest in HVC, since no additional information on any of the variables had value to the decision-maker.

correlations of uncertain variables with the outcome variable, while red bars indicate negative correlations.

## Impact of the Intervention on Biodiversity Conservation Efforts

Biodiversity conservation is linked to income from beekeeping, because conserving the surrounding forests provides flowers for bees to forage on, thus ensuring high honey productivity throughout the year. However, raising incomes for communities within conservation areas, which is the goal of the HVC intervention, does not necessarily mean that these communities will conserve biodiversity [21]. Therefore, for IGAD-BMP to maximize its impact in terms of biodiversity conservation, the program must ensure that rural communities not only generate sustainable incomes from honey production, but also perceive relationships between these incomes and biodiversity conservation, i.e., that high honey production correlates with success in preserving biodiversity. A study on how increased incomes from honey value chain interventions would affect biodiversity conservation outcomes could produce valuable information on the prospective impacts of the intervention on biodiversity conservation.

## Recommendations for Project Monitoring

Although the intervention looks promising, a few areas need to be monitored to reduce the chance of negative outcomes. In Witu, the price of honey has to be monitored to ensure farm gate prices can provide sufficient income to bee-keeping farmers. To reduce the chance of negative outcomes as a result of low farm gate prices, there is a need to link farmers to formal markets, promote collective marketing and invest in value addition. The amount of honey produced per hive should also be monitored, since it is critical in determining whether farmers investing in HVCs are likely to reap net benefits. Adequate farmer training and mitigation of risks, such as theft, land conflicts and fire outbreaks, will help reduce the chance of farmers abandoning honey production. Further, there is a need to identify the minimum number of hives per farmer that would provide substantial increase in income. The implementing agency would then focus on ensuring that most of the farmers have this minimum number to motivate them to continue contributing to HVCs.

In Awer, not much can be done in terms of project monitoring to reduce negative outcomes, since the chance of such negative outcomes was very low. However, the implementing agency can focus on monitoring how these positive outcomes affect biodiversity conservation efforts. There are no clear linkages between increase in income and biodiversity conservation. Clarifying and generating awareness on these linkages, e.g., through research and education on the need to conserve biodiversity to increase honey production, will help the implementing agency to better understand—as well as enhance the impact of—the project on biodiversity conservation.

#### Future Application of the SIE Approach

The costs, benefits and risks of biodiversity-based value chains, and agricultural interventions in general, typically have high levels of uncertainty, especially when considering the longterm and off-site effects of proposed interventions. In light of these uncertainties, the single most important task facing agricultural intervention planners is perhaps to determine the best way to make decisions [22]. Uncertainties and risks have to be quantified, policy preferences clarified and priority measurements for supporting their decisions identified.

The SIE approach can provide plausible solutions to this. The first critical step in the SIE process is to clearly define the decision to be taken. This step often does not receive sufficient attention [23], but it is instrumental in understanding what has to be measured. In this analysis, we spent about 30% of the entire effort on clarifying questions such as whether the implementing agency would cater for the full project implementation cost or whether these costs would be shared with farmers, who would be the actual beneficiaries of the project and whether the implementer was interested in return on investment. Clarifying these questions early in the implementation planning process enabled IGAD-BMP stakeholders to easily identify relevant benefits, costs, risks and external factors that would affect project performance.

The SIE approach allowed identification and prioritization of critical uncertainties in the project. Variables prioritized for further measurement in one of the regions (Witu) appeared easy to quantify with small efforts in data collection. In the other region (Awer), no measurements on any of the variables would change recommendations to farmers and the implementing agency to invest in HVCs. This demonstrates that the SIE approach has potential to substantially enhance the cost effectiveness of decision-supporting measurement campaigns for agricultural intervention planners. The approach presents a clear business case for or against investment in a project where outcomes are uncertain. This can provide critical information to investors aiming to support agricultural interventions with limited resources.

## CONCLUSIONS

Investment in HVC as an economic incentive to conserve biodiversity requires thorough analysis of investment options to maximize returns to different stakeholders. This can be achieved through probabilistic stakeholder-disaggregated outcome projections for biodiversity-based value chain interventions. Clarifying the decision question turned out to be critical for the decision making process. This greatly facilitated the development of a decision's impact pathway, which could then be easily converted into a quantitative decision outcome projection model.

Tailoring interventions to meet economic, social and environmental requirements of rural communities is very important, so decision-makers need approaches that allow holistic ex-ante analysis of investment options. The cost-benefit and risk analyses for biodiversity-based value chains should consider all factors that are relevant for the implementation of the decision, even those that initially appear "intangible" or for which no data are available. Incorporation of local and expert knowledge into decision making using the SIE framework can significantly improve the quality of decisions.

Value of Information analysis can provide indications of what needs to be measured to support intervention decisions. While many uncertainties usually exist in all decisions that affect complex systems, only those uncertainties that are of value to the decision maker should be prioritized for further measurement. This can substantially reduce the cost of data collection aimed at informing decisions.

## STUDY ETHICAL APPROVAL AND CONSENT

The analysis relied on information obtained from literature review, key informants and subject matter experts. This did not require ethical approval as per institutional and national guidelines, but a written and informed consent was obtained from key informants and subject matter experts. The baseline survey outcome in the study background was generated from an open access dataset, where the data was properly anonymized [24]. Written and informed consent was also obtained from all survey participants at the time of original data collection.

## AUTHOR CONTRIBUTIONS

JW: The first author. YK: Contributed to the decision model development. YT: Contributed to the decision model development. GM: Contributed to the decision model development. CM: Contributed to the decision model development. GK: Contributed to the designing of the questionnaire for the baseline survey. Also provided insights on implementation of biodiversity based value chains. JD: Provided expertise on impact of biodiversity based value chains. Also contributed to the editing of the manuscript before submission. JN: Provided vital insights in the framing of the investment decision. Also contributed to the variable estimation process. KS: Provided expertise in decision modeling. EL: Provided expertise in decision modeling and reviewed the manuscript before submission. Also provided insights on the future application of the SIE approach.

#### ACKNOWLEDGMENTS

We acknowledge funding by IGAD's Biodiversity Management Program (IGAD-BMP) and the continuous support of the CGIAR Research Program on Water Land and Ecosystems (WLE). We also appreciate the support accorded to us by the

#### REFERENCES


International Centre of Insect Physiology and Ecology (icipe) through the African Reference Laboratory for Bee Health, led by Suresh Raina. We acknowledge the input by established beekeepers in Mwingi, who provided valuable insights to the modeling team. We also appreciate local communities in Witu and Awer for taking their time to answer questions during the baseline survey.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams. 2018.00006/full#supplementary-material


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Wafula, Karimjee, Tamba, Malava, Muchiri, Koech, De Leeuw, Nyongesa, Shepherd and Luedeling. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Amnesiac Lookback Option: Selectively Monitored Lookback Options and Cryptocurrencies

Ho-Chun Herbert Chang\* and Kevin Li

Department of Mathematics, Dartmouth College, Hanover, NH, United States

#### Edited by:

Francisco Javier Almaguer, Universidad Autónoma de Nuevo León, Mexico

#### Reviewed by:

Daniele Marazzina, Politecnico di Milano, Italy Aurelio F. Bariviera, Universidad Rovira i Virgili, Spain Oleg Kudryavtsev, Russian Customs Academy Rostov Branch, Russia

\*Correspondence:

Ho-Chun Herbert Chang h.herbert.chang.18@dartmouth.edu

#### Specialty section:

This article was submitted to Mathematical Finance, a section of the journal Frontiers in Applied Mathematics and Statistics

> Received: 17 December 2017 Accepted: 19 April 2018 Published: 22 May 2018

#### Citation:

Chang H-CH and Li K (2018) The Amnesiac Lookback Option: Selectively Monitored Lookback Options and Cryptocurrencies. Front. Appl. Math. Stat. 4:10. doi: 10.3389/fams.2018.00010 This study proposes a strategy to make the lookback option cheaper and more practical, and suggests the use of its properties to reduce risk exposure in cryptocurrency markets through blockchain enforced smart contracts and correct for informational inefficiencies surrounding prices and volatility. This paper generalizes partial, discretely-monitored lookback options that dilute premiums by selecting a subset of specified periods to determine payoff, which we call amnesiac lookback options. Prior literature on discretely-monitored lookback options considers the number of periods and assumes equidistant lookback periods in pricing partial lookback options. This study by contrast considers random sampling of lookback periods and compares resulting payoff of the call, put and spread options under floating and fixed strikes. Amnesiac lookbacks are priced with Monte Carlo simulations of Gaussian random walks under equidistant and random periods. Results are compared to analytic and binomial pricing models for the same derivatives. Simulations show diminishing marginal increases to the fair price as the number of selected periods is increased. The returns correspond to a Hill curve whose parameters are set by interest rate and volatility. We demonstrate over-pricing under equidistant monitoring assumptions with error increasing as the lookback periods decrease. An example of a direct implication for event trading is when shock is forecasted but its timing uncertain, equidistant sampling produces a lower error on the true maximum than random choice. We conclude that the instrument provides an ideal space for investors to balance their risk, and as a prime candidate to hedge extreme volatility. We discuss the application of the amnesiac lookback option and path-dependent options to cryptocurrencies and blockchain commodities in the context of smart contracts.

Keywords: options pricing, lookback options, path-dependent options, Monte-Carlo methods, cryptocurrency, smart contracts

## 1. INTRODUCTION

Lookback options have been a prototypical example of exotic options within the financial literature [1]. The option gives the holder the right to buy or sell an underlying asset at any price attained within a specified "lookback" period. The payoff of a lookback call (put) is therefore the difference between the underlying price at maturity and the maximum (minimum) price attained. The trader is thus able to capitalize on the underlying asset as if he sold it at the optimal time. By utilizing only the highest value of the underlying asset in determining payoff, lookback options capture the best case scenarios that people would like to sell at, but often miss due to uncertainty. By corresponding payoffs to the extreme movements of underlying asset prices, lookback options allow for investors to hedge against or invest in volatility.

Goldman et al. [2] [3] are one of the first who mentioned the lookback option within the financial literature. Their paper examines the properties of path-dependent European options under Black-Scholes assumptions. They laid out three primary motivations: to minimize regret, to live out the "fantasy" of buying low and selling high, and lastly, to use knowledge of the range to expand an investors opportunity set. This last condition is only applicable within the realistic market setting rather than a frictionless one. The nature of lookback options also lets investors guard against the behavioral flaws in people. The nature of lookback options also lets investors guard against the behavioral flaws in people. Because the payoffs of lookbacks are path dependent in a way to capture the effects of the best prices, lookbacks can hedge not only quantitative variations in the market, but irrational regret of the human mind.

Lookback options benefit the holder given greater volatility in the market. Thus, it is an effective instrument for hedging against large price movements and reducing the risk of destabilizing events that may cause markets to either rise or fall. These include political election outcomes or other unforeseen chance occurrences such as market anomalies [4]. Also among its potential uses is to hedge exposure in cryptocurrency markets, where extreme volatility has increasingly become a menace as investors look for ways to mitigate such risk.

However, the lookback option's strong reduction to risk exposure requires a sufficiently high premium, which reduces the option's demand. The proposition of this paper recognizes that one may not have to look back upon the entire associated time period; merely looking back upon a portion suffices to capture most of the max potential payoff of a full-time lookback option, most of the time. This conclusion allows for investors to adjust their risk using lookbacks of differing lookback period lengths by purchasing shorter time lookback options at a discount in order to increase return on investment. This strategy suggests the possibility of event-based trading whereby lookback periods may be chosen to correspond to market shaping events in hopes of extraordinary profits.

As the potential payoffs and insurance capabilities of an option increase, so does its premium. Thus, there is a desire for an instrument that limits risk exposure, but at a lower price. Here is where the proposed "Amnesiac Lookback Option" may find a tradeable niche, a variant of the conventional lookback option whose lookback time periods are restricted to a chosen subset of periods within a specified time interval at the time the contract is created. If premiums for lookbacks may be reduced, it could also open the door for their widespread use, particularly in high volatility cryptocurrency markets. The purpose of this study is therefore 2-fold: to propose a way to make the lookback option cheaper and more practical, and suggest the use of its properties to reduce risk exposure in cryptocurrency markets and correct for informational inefficiencies surrounding prices and volatility. Additionally, we are interested in the following research questions. What periods of the amnesiac lookback option should be selected to maximize its final payoff? What is the probability that the amnesiac lookback option may produce greater profit than the standard lookback option? Lastly, would this option have practical applications in the cryptocurrency market?

The rest of the paper is organized in the following way: The rest of the introduction provides more context on the varieties of lookback options and applications to cryptocurrencies. Section 1 continues to discuss types of lookbacks and its application to the cryptocurrency market. Section 2 discusses existing literature on partial lookback options, lookback options with diluted premiums, and past efforts to price them. It then discusses our methods of simulation. Section 3 presents the results of three types of numerical simulations. First, a demonstration of the general properties of floating amnesiac lookback options for calls, puts and spread options are shown. Second, the payoff of fixed strike amnesiac lookback options is shown and compared to the fixed strike. Third, a study of equidistant monitoring vs. two forms of random monitoring are presented. Section 4 prices different cryptocurrency using the algorithm in conjunction with smart contracts.

## 1.1. Fixed vs. Floating Lookback Options

Lookbacks exist in many different varieties but can be classified into two broad categories: fixed and floating strikes. The strike price of fixed strikes is indicated in the initial contract, while for floating strike lookbacks the strike price is the minimum (call) or maximum (put) of the underlying asset. Additionally, there are spread lookback options with payoffs equal to the difference between the maximum and minimum prices attained by the underlying. The payoffs are summarized in **Table 1** where M<sup>T</sup> <sup>0</sup> = max{Si} T i=0 and m<sup>T</sup> <sup>0</sup> = min{Si} T i=0 . T denotes the time of maturity.

Suppose an investor decides to take the long position for the next 2 months. However, the price of the stock drops unexpectedly within the last 4 days of closing by 10%. Not selling the stock early becomes a reason of regret. Similarly, the stock may drop in the span of 2 days before rising quickly, making the investor regret not buying the stock a little later. Thus the floating strike call lookback is useful for market exit and the fixed strike call for market entry.

## 1.2. The Partial and the Amnesiac

Due to the high premiums of lookback options and perhaps other factors, lookback options are not sold at particularly high volumes within OTC trades [5]. It has been suggested that partial lookback options may allow lookbacks to be more tradeable. One

TABLE 1 | Payoffs for floating and fixed strike lookback options.


such way is to introduce a factor λ that decreases the effect or payoff of the lookback, proposed by Conze and Viswanathan [6]. For instance, the payoff of such an instrument is:

$$\lambda \left( \max \{ \mathcal{S}\_i \} \_1^N - \mathcal{S}\_T \right) \tag{1}$$

If λ equals 1 then payoff is identical to a standard lookback option. Another form of the partial lookback, referred to as the fractional lookback, restricts the lookback periods to a continuous subset of days. Finally there is also the window lookback option, where separated continuous periods are monitored [7].

The amnesiac lookback option generalizes partial and discrete lookback options that are not linearly reliant on the payoff of the standard lookback option and is named such in response to Heynen and Kat's paper Selective Memory. The amnesiac lookback option is a lookback option whose lookback periods belong to a pre-specified subset of periods from a given time interval, determined at the moment of contract creation.

The amnesiac lookback option is thus defined as follows. Given a full, discrete lookback option with N total lookback periods, an amnesiac option is defined by a subset of periods of the full, discrete lookback option, denoted as A. The payoff scheme is the same as the standard lookback, while the extreme values of the option are restricted to A. This means, every full, discrete lookback option is in fact an amnesiac option of a full lookback with greater periods. This extends previous definitions of partial lookbacks to beyond equidistant intervals.

#### 1.3. Application to Cryptocurrencies

Lookback options are extremely suited to hedge against volatility in general, whether the underlying asset surges in value, or whether the underlying asset declines in value. Unpredictability is therefore the prime feature of an asset that would drive demand and usage for a lookback option based upon that asset. In 1982, the Mocatta Metals Corporation issued one of the first "lookbacks," that allowed a trader to buy gold at the lowest price attained within a period. In the context of modern day financial markets, it would seem lookbacks could have high potential in hedging investments involving the popular digital gold cryptocurrencies. A pioneering innovation within currency markets, instruments such as Bitcoin and Ethereum may represent the future of monetary exchange. Given its relative technology security, explosion in value, and increasing acceptance of legitimacy, Carrick [8] has even suggested that cryptocurrencies like Bitcoin could be used as a complement to fiat currencies in emerging markets. At the same time, high volatility has become the premier feature of cryptocurrency markets, which has made investment risky. A Bitcoin was worth \$2 in 2011—and exploded to \$4,000 in 2017 [9]. During that time, volatility for daily returns would regularly exceed 10%, even skyrocketing to 16% [10]. Bitcoin once lost 75% of its value over 2 years, and then rising 2,000% within the next two.

Past literature has found evidence of time series characteristics and long memory behavior in Bitcoin markets, both in regards to pricing and in return volatility. Bariviera et al. [11] explores some of these inefficiencies by calculating Hurst exponents via the Detrended Fluctuation Analysis method for certain time windows of Bitcoin return and return volatility data. Bariviera finds evidence that from 2011 to 2014, daily Bitcoin returns had long-term positive autocorrelation with previous returns and that return volatility had long memory throughout the entire time period of 2011–2017 [11]. Bariviera et al. [11] utilize this same methodology to further show that intraday returns before 2014 exhibit long range memory as well. Additionally, Urquhart [12] finds that Bitcoin prices exhibit odd behaviors through the entire time period May 2012 to April 2017, with over 10 percent of prices ending with decimal digits of 00. One, two, three, five, and ten days before a round number from rising prices, returns are positive and statistically significant. Meanwhile, prices after a round number show no such predictable behaviors. Given that the lookback option utilizes the time series behavior of an asset's prices and volatility to determine its payoff, it is a suitable candidate for rectifying such market inefficiencies. The amnesiac lookback option is particularly appropriate because lookback period selection allows for targeted correction of unusual events, such as those found in Urquhart [12], and can more easily eliminate arbitrage opportunities and bring markets into more efficient states.

Ultimately, Bitcoin aims to be a type of fiat money; it has no value backed with consumable goods and its value comes from the minds of its investors and from the financial environment in which it occupies. Increases have been driven by hopes of future value, or in other terms, heavy speculation. Unexpected events therefore have rippling effects on the cryptocurrency market and create a volatility with few equal comparisons. China's decision to cease its bitcoin exchange in September of 2017 sent Bitcoin into a downward spiral, and as governments and regulators venture forth into the frontier that is the cryptocurrency market, events that will shock the market are sure to take place [13]. Even so, investors are increasingly accepting of cryptocurrencies as attractive investment propositions outside of speculation. Bouri et al. [14] find that Bitcoin can act as a hedge against market uncertainty in situations, specifically in short time horizons under extreme bear or bull market regimes, or when uncertainty is either very low or very high. Bouri et al. [15] overall minimizes the usefulness of Bitcoin as a general hedge, but still finds evidence that investors put money into Bitcoin when Asian stocks experience extreme down movements.

In July 2017, the U.S. Commodity Futures Trading Commission granted federal approval to LedgerX LLC for the ability to trade and exchange options based on cryptocurrency [16]. This happening marked the birth of the first federally regulated platforms of cryptocurrency options trading, and moved items such as Bitcoin, Ethereum, and even Dogecoin closer into the realm of financial legitimacy. Shortly after on December 1, 2017, the U.S. Commodity Futures Trading Commission announced the offering of self-certified derivatives from three financial firms: bitcoin futures from the Chicago Mercantile Exchange Inc. and the CBOE Futures Exchange and binary options from the Cantor Exchange [17]. In contrast with the effect of regulatory restrictions, this decision has driven up bitcoin prices from \$5,000 prior to the announcement to above \$11,000 by December 05, 2017 [18]. The issuance of new financial instruments can thus feedback into the price of its underlying.

More importantly, the integration of smart contracts into blockchain technology has vastly expanded the possibilities of options trading, allowing traders great flexibility in designing their own options [19–21]. Smart contracts are computer protocol that allow for trade and exchange without the need for an intermediary, and because writers code their own contracts, blockchain technology can facilitate the trade of highly customized contracts often seen on the OTC market. Both Bitcoin and Ethereum, for instance, support programming languages that allow for the creation of custom smart contracts [22, 23]. Amnesiac lookback options could conceivably exist within this structure, and be traded as a smart contract.

Recent events lent credence to the prospect of a prolific smart contract exchange. In March 2018. The state of Tennessee signed into law a bill that recognizes smart contracts as having legal power [24], providing a pathway for other states to follow and further securing the legitimacy of these digital arrangements. Blockchain based smart contract firm Hedgy has also created irrefutable and unalterable "Smart Futures" that can enforce digital obligations and streamline settlements [25]. As more and more investors begin creating and trading new smart contracts, blockchain may 1 day even "democratize" the OTC market and open a plethora of new smart instruments to be exchanged.

Lookback options therefore emerge as an attractive and easily implementable product for currency speculators. As price fluctuation is the determining factor of payoff, such lookback instruments would be highly priced in the cryptocurrency market and be a key counterweight to the speculative risk of investors. Simple to understand and able to exist as smart contracts, lookback options and other path-dependent options could serve as a unique tool for traders with restricted access to dynamic trading strategies, and in the process introduce greater flows of capital from these sources to promising investments.

## 2. MATERIALS AND METHODS

#### 2.1. Past Pricing Model for Lookbacks

Research with lookback options fall within pricing methods for path dependent exotic options, such as the barrier option which also rely on extreme value processes, and the Asian option which relies on the average price of the underlying. Goldman introduced the instrument into the financial literature in 1979 [2, 3]. Following Black, Scholes and Merton's pricing of vanilla options, Black-Scholes assumptions have been extended to price exotic options. This is defined by at least one asset with price S that moves under Geometric Brownian Motion with constant drift and volatility. The underlying price follows:

$$dS = \mathcal{S}\mu dt + \sigma \mathcal{S} dX \tag{2}$$

where the underlying has drift µ and volatility σ. In the riskneutral measure, the stock price at a given time t and final time T is given below:

$$S\_{t+1} = S\_t e^{\left(r - D - \frac{\sigma^2}{2}\right)\Delta t + \sigma \epsilon \sqrt{\Delta t}} \tag{3}$$

$$S\_T = S\_0 e^{(r - D - \frac{\sigma^2}{2})T + \sigma \epsilon \sqrt{T}} \tag{4}$$

r denotes the risk free rate, σ the volatility, D denotes dividends in the case of stocks, and ǫ ∼ N (0, 1), the normal distribution.

Like the Vanilla European Option, an analytic pricing formula has been shown using martingale methods [26], with the price of a call given as:

$$C\_t = \left\{\Phi(a\_1(\mathbb{S}, m)) - me^{-r\tau}\Phi(a\_2(\mathbb{S}, m)) - \frac{\mathbb{S}\sigma^2}{2r} \{\Phi(-a\_1(\mathbb{S}, m))\} \right. \tag{5}$$

$$-e^{-r\tau}(\frac{m}{\mathbb{S}})^{\frac{2r}{\sigma^2}}\Phi(-a\_3(\mathbb{S}, m)) \} \tag{5}$$

$$P\_l = -\mathbb{S}\Phi(-a\_1\{\mathbb{S},M\}) + Me^{-r\tau}\Phi(a\_2\{\mathbb{S},M\}) + \frac{\mathbb{S}\sigma^2}{2r}(\Phi(a\_1\{\mathbb{S},M\})) $$

$$-e^{-r\tau}(\frac{M}{S})^{\frac{2r}{\sigma^2}}\Phi(a\_3\{\mathbb{S},m\})) \tag{6}$$

where M denotes the running maximum at time t, m denotes the running minimum,τ = T − t with T the time of maturity, and 8 the standard normal cumulative distribution function given as

$$\Phi(\alpha) = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\alpha} e^{-\frac{x^2}{2}} dx$$

With L being a dummy variable, the variables denote:

$$\begin{aligned} a\_1(\mathbb{S}, L) &= \frac{\ln \frac{\mathbb{S}}{L} + (r + \frac{1}{2}\sigma^2)\tau}{\alpha \sqrt{\tau}}\\ a\_2(\mathbb{S}, L) &= a\_1(\mathbb{S}, L) - \sigma \sqrt{\tau} \\ a\_3(\mathbb{S}, L) &= a\_1(\mathbb{S}, L) - \frac{2r\sqrt{\tau}}{\sigma} \qquad \text{for } L > 0, S > 0 \end{aligned}$$

For fixed strikes, Conze and Vishwanathan [6] used known properties of maxima and minima distributions to price the lookback call, under the continuous case of lookback periods.

$$\mathbf{C}\_{t} = \begin{cases} \mathbf{S}\_{0}\Phi(b(T)) & -e^{rT}K(b(T) - \sigma\sqrt{T} + e^{-rT}\frac{\sigma^{2}}{2r}\mathbf{S}\_{0}\big(e^{rT}\Phi(b(T))\big) \\ & -\left(\frac{\mathbf{S}\_{0}}{K}\right)^{\frac{-2r}{\sigma^{2}}}\Phi(b(T) - \frac{2r}{\sigma}\sqrt{T})\big) & K \ge M \\ e^{rT}\langle M\_{T}^{0} - \qquad K\rangle + \mathbf{S}\_{0}\Phi(b'(T)) - e^{-rT}M\_{T}^{0}\Phi(b'(T) - \sigma\sqrt{T}) \\ & + e^{-rT}\frac{\sigma^{2}}{2r}\mathbf{S}\_{0}\big(e^{rT}\Phi(d\_{T}^{\prime}) - \left(\frac{\mathbf{S}\_{0}}{M\_{T}^{0}}\right)^{\frac{-2r}{\sigma^{2}}}\Phi\big(b'(T)\big) \\ & -\left(\frac{2r}{\sigma}\sqrt{T}\right) & K < M \end{cases}$$

where

$$b(T) = \frac{\ln\frac{S\_0}{K} + (r - D + \frac{\sigma^2}{2})T}{\sigma\sqrt{T}}$$

$$b'(T) = \frac{\ln\frac{S\_0}{M\_0^T} + (r - D + \frac{\sigma^2}{2})T}{\sigma\sqrt{T}}\tag{7}$$

The similar case for the put is derived in their paper, and resembles the form of the float-strike option in Equation (5). These analytic solutions assume continuous monitoring and updates to the price of the underlying, while in reality most options are monitored at discrete intervals. Thus, the analytic solutions are typically overpriced as compared to the prices of discretely monitored lookbacks, and has spurred two areas of study: how to either price discrete options more accurately, and how to correct errors of continuous monitoring.

Historically, notable papers that have priced discrete lookback options include those by Heynen and Kat using the Black-Scholes model [27], Babbs using a binomial valuation [28], and Cheuk and Vorst [29] using a binomial model. More recent pricing paradigms utilize stochastic volatility [30] and the Lévy Model [31] (more citations here). Cheuk and Vorst in particular consider the effects of observation frequency on price using equal lookback intervals.

The binomial options pricing model first described in Cox et al. [32] offers a conceptual method with which the price of an amnesiac lookback options may be determined. Under the assumptions of the binomial tree model, the price of an underlying stock S<sup>t</sup> at discrete time t may strictly exhibit one of two behaviors. It may either increase in value to Stu with probability p, or decrease in value to Std with probability q where u = e σ √ t , d = 1/u, p + q =1, and sigma is a measure of the stock's volatility. The risk-neutral probability p is determined to equal e (r1t)−d u−d where r an exogenously determined risk free return. The lookback option can then be priced recursively using a lattice.

This pricing model was then extended to be more efficient by Conze and Vishwanathan [6], Babbs [28], and Cheuk and Vorst [29]. The binomial pricing method lends itself more easily to pricing at discrete intervals. With too few time intervals, the approximation for the standard lookback option becomes inaccurate. Yet as the number of periods in the binomial tree increases, so does the amount of computations. Past work involves increasing the efficiency of these computations using combinatorial properties of the binomial tree. Babbs [28] uses a change of numeraire approach, similar to Hull and White [33], to price both new and existing lookback options. Cheuk and Vorst formulate a modified tree, where the computation only relies on time and the difference between u and d. Let this variable be called j such that:

$$\mathcal{S}\_t = \left(\min\_{t \le T} \mathcal{S}\_t\right) \boldsymbol{\mu}\_n^j \qquad \mathcal{S}\_t = \left(\max\_{t \le T} \mathcal{S}\_t\right) \boldsymbol{\mu}\_n^{-j} \tag{8}$$

for the call and the put respectively. Then the price for a floating strike lookback can be given as C fl <sup>n</sup> = S0Vn(0, 0), where:

$$V\_n(0,0) = \sum\_{j=0}^n (1 - \mu\_n^{-j}) P(V\_{0,n}^\* = 0, V\_n^\*, n = j) \tag{9}$$

The notation star such as in V ∗ m,n denotes for n fixed and 0 < m < n. The probabilities are given by:

$$P(V\_{0,n}^\* = 0, V\_n^\*, n = j) = \sum\_{k=j}^l \Lambda j, k, m(1 - q\_n)^k q\_n^{m-k}$$

$$\text{for} \qquad \Lambda j, k, m = \binom{m}{k-j} - \binom{m}{k-j-1} \tag{10}$$

if k > j and 3j, k, m = 1 is k = j.

While both the modified tree model and the analytic formula give continuous valuation results, and it is well known the binomial tree model converges to the continuous price. the pricing formulae is still related to discrete lookback pricing. Cheuk and Vorst extended the total number of lookback periods, then selected an equal-distant subset to model the discrete monitoring periods. Similarly, the amnesiac lookback option is defined on a subset of monitoring periods.

Upon noting the discrepancy between continuous and discrete monitoring pricing assumptions, researchers continued to improve methods of pricing discrete options. Aitsahlia and Le Lai [34] uses the random walk duality to simplify recursive integration. Petrella and Kou [35] uses the Laplace transform to make pricing more efficient, and Broadie et al. [36] uses approximation adjustments to the continuous case to fairly price the option's discrete counterpart.

Contemporary researchers derive pricing under the framework of Lévy Processes, and as a natural extension, the Wiener-Hopf factorization is implemented as it gives the distribution of functionals within a random walk. Fusai et al. [37] combined this with Spitzer's identity for a general method of pricing discretely monitored exotic options, with an explicit formula for the fixed strike option. However, the method assumes equidistant monitoring windows. Boyarchenko and Levendorskii [38] similarly implement the Wiener-Hopf factorization but correct for the discrete monitoring dates using a Laplace transformation. Hieber [39] presents the general pricing of exotic options using the Fourier Transform and under Black-Scholes assumptions and regime switching model. This is not just useful for lookback options but for digital options and barrier options. However, the pricing scheme is for continuously monitored options. Feng and Linetsky [40] on the other hand, implements Hilbert transformations sequentially and presents an interesting method of computing maxima and minima, thereby applying it to the valuation of path-dependent options.

#### 2.2. Method of Simulation

The definition of a amnesiac option in section 1 can be formalized as follows. Let N denote the total number of periods of a full lookback option. Let A denote a selected subset. Given a sequence of asset prices{Si}, define M<sup>T</sup> 0 the maximum and m<sup>T</sup> 0 the minimum defined on the subsequence where i ∈ A. The payoff of the amnesiac option is then the same as **Table 1**.

To price the amnesiac lookback option, we utilize a Monte Carlo pricing simulation under the assumptions of the Black-Scholes world, including but not limited to, its modeling parameters, risk-free rate, dividends, and volatility. We assume geometric Brownian Motion under a risk-neutral measure as described in Equation (4):

$$\mathcal{S}\_T = \mathcal{S}\_0 e^{(r - D - \frac{\sigma^2}{2})T + \sigma \epsilon \sqrt{T}}$$

We assume dividends D = 0. We simulate this process using a Gaussian random walk. The algorithm is shown explicitly in Algorithm 1:


Confidence Interval = Fair Price ± (t)(SE)

where t is the corresponding t value for the confidence level. The updating of the asset price follows a Gaussian Random walk in Algorithm 2 where dT is the total time T divided by the number of selected days p:


A denotes how the days are chosen. Our experiment features three forms of period selection: uniform random with fixed endpoints, completely random, and equal distant selection. For the amnesiac option with fixed endpoints, the first and last indices are first chosen, then indices 0 < i < N are chosen using the uniform distribution. For true random, indices are randomly selected under the uniform distribution with no guarantee of the first one being selected. For equal distant period selection, indices are chosen, and rounded down when not divisible. We make an assumption that on average, the population's selection of periods will be uniform random.

The parameters are shown in **Table 2**.

For each time step, the maximum and minimum value is updated. We analyze the pricing of floating and fixed amnesiac lookback options, along with the variants discussed in **Table 1** TABLE 2 | Simulation parameters.


(call, put, and spread). Then, the payoffs under equal interval sampling and random sampling are considered.

A useful function that appears in analysis is the general Hill Function. the Hill Function in the form:

$$V(X) = \frac{V\_{\text{Max}}X^h}{K^h + X^h} \tag{11}$$

The Hill Function has its roots in chemistry, where it is commonly used to model the kinetics of substrate reaction rate, where the x-axis is the substrate concentration and the y-axis the rate of reaction. VMax is the maximal saturation rate, K denotes the value of x required to attain half of VMax, x is the density of a substance, and h is the hill coefficient, as an exponential of K and h.

Analogously, increasing lookback periods yield diminishing returns to payoff and the components of the Hill Function give intuition to the process itself. Firstly, Vmax is theoretically the maximum payoff, or the price of the standard discrete lookback option. Secondly, as we increase the starting stock price, both Vmax and C scale linearly. On the other hand, both K and h remain constant and depend on interest rate r and volatility σ. The Hill Curve can thus price amnesiac lookbacks independent of the initial asset price, keeping interest and volatility constant. If the price of the standard lookback option Vmax is known, and interest rate and volatility incorporated into K and h, then pricing can be done very efficiently using the Hill Equation. Additionally, the curve can be adopted to any starting stock price, since the initial stock price only influences Vmax. The ratio between the √

price evolution and initial price is simply <sup>S</sup>0<sup>e</sup>

$$\frac{S\_0 e^{(r-d-\frac{\sigma^2}{2})T + \sigma\sqrt{T}}}{S\_0} = $$

e (r−d− σ 2 2 )T+σ √ T , hence when pricing under the same interest and volatility the payoff can simply be scaled.

The total number of simulations per parameter set is 2,000,000, the fair price being the arithmetic mean of the payoffs, with a target maximal standard error of 0.001.

## 3. RESULTS

We denote C<sup>n</sup> and P<sup>n</sup> as the amnesiac option with n chosen monitoring periods. When unspecified, the monitoring regime is assumed to be random period monitoring with fixed endpoints.

## 3.1. Floating Strike Amnesiac Option

To reiterate the definition of an amnesiac lookback option, its payoff scheme is as described in **Table 1**, where only the maxima and minima of a subset of periods are recorded. We use T and t to

denote the continuous case and N and n be the discrete case. Let N be the total number of periods and A the ordered set of indices:

$$A = \{n\_a, n\_b, n\_N\} \qquad \forall A \subset N$$

The amnesiac is defined on any k selected monitoring periods.

**Figures 1**, **2** show the payoff vs. the number of days selected. As the number of periods look backed on increases, the payoff approaches that of the standard lookback option. Yet it does not require many selected periods to capture most of the payoff of the full 100 selected days. This can be attributed to the diminishing marginal effect of additional selected periods; each period selected after the first contributes less additional payoff than the previous. For each data point, we took the mean of 2, 000, 000 simulations, and given a calculated Standard Deviation of around 14.54, this gives an error of around 0.01 in absolute terms.

The sample errors are summarized in **Table 3**.

**Figure 1** fits the Hill function as previously described, and **Figure 2** extends this to include volatility. For this particular case, TABLE 3 | Average of sample errors from 2,000,000 simulations.


the formula is observed to be:

$$V(X) = \frac{30.96X^{0.62}}{1.18^{0.62} + X^{0.62}}$$

Kernel density estimators further tell the differences of portfolios constructed from different amnesiac. **Figure 3** shows the distribution of the call payoffs. The peak clustered at zero denotes the likely-hood of the option not being exercised, which yields a payoff of zero. Comparisons between instruments vary case-bycase, so it is useful in comparing expected returns in a pairwise fashion.

**Figure 4** shows the expected returns for the number of monitoring periods n equal to 99, 7, and 2. We define returns by the payoff minus the price of the option, which is the mean. Let x<sup>i</sup> denote a single payoff from the simulation, then let the price of the call be defined as the average X¯ then the returns are defined by:

$$\mathcal{R} = \varkappa\_i - \bar{X}$$

The payoff distribution of n = 7 for R > 0 is approximately the same as the distribution of n = 99. In fact, they converge in the negative region, hence their probability of R > 0 is approximately equal. Next, note that the downside of n = 99 is greater than that of n = 7. This means that the downside for C<sup>7</sup> is more limited than C99, even if their conditional expectation for the downside is equal or even worse. This indicates that the amnesiac lookback, and partial lookbacks with less monitoring periods in general, can be used to limit the downside effects.

#### 3.1.1. Three Types of Payoffs

**Figure 5** shows the three payoff types as described in **Table 1** for floating strikes, denoted in the first column. The sum of the fairprice of the call and put is equal to the spread option, which is why it is often described as a "double lookback option."

**Figure 5** also shows that the average maximal and minimal value is independent of the average final value. This observation can be illustrated as follows, given the payoffs in **Table 1** and **E** the expectation operator.

$$\begin{aligned} \mathbb{E}(call) + \mathbb{E}(put) &= \mathbb{E}(\mathbb{S}\_{max} - \mathbb{S}\_{T}) + \mathbb{E}(\mathbb{S}\_{T} - \mathbb{S}\_{min}) \\ &= \mathbb{E}(\mathbb{S}\_{max} - \mathbb{S}\_{T} + \mathbb{S}\_{T} - \mathbb{S}\_{min}) \\ &= \mathbb{E}(\mathbb{S}\_{max} - \mathbb{S}\_{min}) = \mathbb{E}(spread) \end{aligned} \tag{12}$$

These results show remarkable resemblance to the rate of convergence for exotic prices described in Dai et al. [41]. Indeed the prices produced by varying the number of total periods

N follows a similar curve, and produces an upper bound for the pricing. Furthermore, since the distribution of stocks is lognormal, if the interest rate r = 0, then the put will be of greater value than the call because the maximum will fluctuate absolutely more than the minimum. To see this suppose the simple case in the binomial model. Then for fixed n, Se<sup>n</sup> √ <sup>σ</sup> − S is greater than S − Se−<sup>n</sup> √ σ . However, due to the positive risk-free rate the price of the call is greater as shown in **Figure 5**.

For illustration, consider the kernel densities of each of the three instruments shown in **Figure 6**. First, note the distribution of the spread option is lognormal, shown in the blue curve. In other words, the distribution of the range of a geometric Brownian motion is lognormal, due to the range of brownian motion being distributed normally.

Characteristic of the lognormal distribution is its fatter left tail. This corresponds to the shape of the put option, minus the peak created by the cluster of 0s. Additionally, the tail of the put is platykurtic in comparison to the tail of the call. In sum, once the 0 payoffs are eliminated, the distribution of the put corresponds to the left tail of the spread, and the distribution of the call corresponds to the inverted right tail.

#### 3.2. Fixed Strike Options

The analysis of fixed strike amnesiac options first requires an analysis of how in-the-moneyness affects option payoff. In **Figure 7** the call price decreases as the strike price increases, as the more out-of-the-money the option gets, the more unlikely the maximum of the underlying will be greater than the strike. The

put option in **Figure 8** demonstrates the opposite relationship, as it relates to the minimum distribution. The payoff diagrams for the call show a linear relationship while in-the-money, then decaying exponentially afterwards. The put shows the inverse. Since the maximum (minimum) distribution for given number of periods is fixed, as the strike price decreases (increases) the difference is linear. This linear relationship is preserved as the number of selected periods in the amnesiac option are reduced, as seen in the diagrams.

When the number of monitoring periods is equal to 0, then the maximum is simply the last stock price at time N. Thus, the payoff is the final stock price minus the strike.The shape of this curve for selecting n periods is therefore a simple, vertical translation. Through these observations, we have demonstrated the bounds for the amnesiac option as shown in **Figure 9**. It is bounded by the standard lookback from above, and the vanilla European options from below.

## 3.3. Random vs. Equal Monitoring

One of the most important questions this paper addresses is if the choice of period matters in the payoff. What we find is randomly choosing lookback periods yields a worse payoff than equally spacing the days to look back on, as shown in **Figure 10**. The differences are shown in **Figures 11**, **12**.

There are two ways to monitoring randomly. We defined N equidistant points between time 0 and T, making it discrete. Suppose we are monitoring k points from the Gaussian Random walk. Under true random monitoring we would select k random points from the set of periods {0, 1, 2, ..., N − 1}. In contrast, the fixed-end random monitoring selects the first period n = 0 then randomly chooses k − 1 periods from indices 1 through N − 1. Since the price of an option increases with the duration of the contract, fixing the first and last period ensures an equal time length for comparison. For the true random monitoring, the time to maturity would be the minimum of index.

For illustration, suppose k = 2. For fixed-end random sampling the first period is chosen while the second can be any point between 1 and N − 1. More simply, both endpoints are fixed, and we are choosing a point between the beginning and end. Numerical results in **Figure 13** on average selection of the midpoint as more optimal. True random sampling produces weak payoffs due to shortened maturities. However, even when endpoints are always chosen to lookback on, the payoff for equidistant sampling is superior. Although this suggests the optimal strategy is to equally space the selected periods, this does not take into account sudden volatility swings from extreme, singular events.

The intuition for this toy experiment can be explained as follows. Since both endpoints are fixed, we are only interested in choosing the point in-between. If the point chosen is less than or equal to S<sup>0</sup> or SN, then the payoff is the same as the case of choosing one point, which is S0. For the third point S<sup>i</sup> to be of value it must satisfy S<sup>i</sup> > S<sup>0</sup> and S<sup>i</sup> > SN. Suppose n is very close to the first index 0, then S<sup>n</sup> = S0e ǫ for ǫ that falls in the upper half of the Gaussian distribution defined in Equation (2). The difference between S<sup>0</sup> and S<sup>n</sup> is thus small, and the payoff will be close to that of the case of just choosing the first point. On the other hand, suppose n is close to N. Then S<sup>N</sup> = Sne ǫ . Since the payoff is just the difference between the max and final value, the payoff is diminished.

Simply put, given that the chosen point S<sup>n</sup> is greater than S<sup>0</sup> and SN, then choosing too close to the beginning means not giving enough distance for the maximum to rise, and choosing too close to the end means diminished payoff due to proximity with SN. Intuitively, time value is then the explanatory factor that determines how much marginal payoff a selected period gives, where the value of a point is proportional to the time value associated with the distance between it and the closest adjacent lookback period.

We now demonstrate the above claim mathematically. For simplicity we fix the first period and denote j = k − 1 as the number of lookback periods. Then let N denote the total number of equally spaced periods of the standard lookback with total time T, and let A be the set of selected periods. Define function γ that gives the minimum distance

$$\gamma(i) = \min \{ d(i, k) \quad \forall k \in A, k \neq i \} \tag{13}$$

Then the totals of time value is given as

$$\Gamma(A) = \sum\_{i \in A} \mathcal{V}(i) \tag{14}$$

where the value of Ŵ for equidistant sampling is <sup>j</sup> j+1 . We show this value is maximal.

Consider the base case of n = 1. If i<sup>0</sup> = N 2 the only element in A0, then Ŵ(A0) = N 2 . Compare to this to A<sup>1</sup> whose only element i<sup>1</sup> = N <sup>2</sup> <sup>+</sup> <sup>ǫ</sup> for <sup>ǫ</sup> <sup>∈</sup> <sup>Z</sup> <sup>+</sup>, ǫ < <sup>N</sup> 2 . Since d(0, i1) > N 2 , then d(i1, N) = N <sup>2</sup> <sup>−</sup> ǫ < <sup>N</sup> 2 . Thus Ŵ(A1) = N <sup>2</sup> − ǫ < Ŵ(A0). By symmetry, any i<sup>1</sup> greater or less than <sup>N</sup> 2 is worth less. Cases for larger j can be shown in a similar fashion by assuming the first point greater than <sup>N</sup> j . Thus our claim that lookback period value is proportional to the time value of the next closest period remains true.

What this conclusion implies in real situation trading is that it would be inefficient to cluster lookback points around a presumed volatile event in order to capture a price extrema, even if volatility spikes seem certain, because the marginal value of each period chosen beyond the first in the locality diminishes rapidly due to the decreasing proximity of neighboring points. Prices are less likely to change a significant amount between two points in time the closer the two points in time are.

From the framework of analysis, suppose we sub-divide Geometric Brownian Motion into N subintervals defined by select indices. Then we are interested in the endpoints [i, i + 1] as our monitoring dates. If these subintervals are of equal length, then the probability of the true maximum being captured is equal, since within the subinterval, the maximum is equally likely to fall on any point. Thus, the probability of the maximum being captured by the neighborhood around the endpoints diminish as subintervals get larger. The maximal neighborhoods are when endpoints are equally spaced.

This has been explored in the risk theory and queuing theory literature, which assumes expected maxima calculations on

FIGURE 9 | The payoffs of the upper and lower bounds for the call and put respectively. The area between demonstrate the customizable area of the instrument payoff.

equidistant sampling. For further references consider papers by Janssen and Van Leeuwaarden [42, 43] on equidistant sampling of Brownian Motion and Gaussian Random Walks, and Alfi's work on exacting roughness in finite random walks [44]. Under equidistant sampling, the expected maximum shown using the Spitzer Identity which gives the joint distribution of partial sums and their maximal sums for a collection of random variables:

$$\mathbb{E}\max\_{n=0,\ldots,N} B(\frac{n}{N}) = \sum\_{n=1}^{N} \frac{1}{n} \mathbb{E}\max\{\frac{n}{N}\}\tag{15}$$

Now let the process X<sup>n</sup> be defined as

$$X\_n = \sum\_{i=1}^n \delta\_{xi} \tag{16}$$

After computation [44], the expected maximum value can be expressed as the following

$$\mathbb{E}(M\_k) = \lim\_{s \to 1} \frac{\mathbb{E}(s^{M\_k}) - 1}{\ln(s)} = \sum\_{i=1}^k \frac{\mathbb{E}(|X\_i|)}{2i} \tag{17}$$

Adapting these results to Gaussian random walks with unequal steps may produce approximations for the payoff of lookbacks and other exotic options.

If the events of the world are analogous to random occurrences, then nonrandom selection of lookback periods is akin to event-based trading strategies. Specifically in regards to the amnesiac lookback option, such a strategy leads to lower payoffs. This conclusion implies that given general market efficiency, event-based period selection strategies of lookback periods should be unprofitable.

## 4. DISCUSSION

These results can be immediately applied to pricing cryptocurrency. Under the Black-Scholes model, the implied volatility is usually calculated using the current market price of the option. However, because this is a new option being proposed, and since both futures and vanilla options are still new derivatives on the market, it is not possible at this time to even price them with rudimentary Monte-Carlo simulations.

Instead, we use the historical volatility as a proxy and use the 3-month treasury bond as the risk-free rate. **Tables 4**, **5** show the prices of amnesiac calls and puts based on historical data in a 3 and 2 month window, respectively. These are then shown

as a ratio of the initial price in **Figure 14**. Data is taken from Coinmarketcap.

The last row in **Table 5** shows a Bitcoin amnesiac option priced at the implied volatility of Bitcoin calls in December. Interestingly, it is less than the minimum volatility of Bitcoin itself. When the implied volatility is less than historic volatility, this indicates that market sentiment benefits option buyers. In other words, the market expects lower volatility in the future, and as a result, using historical volatility as a proxy is prone to overpricing.

To verify this is what we expect, we observe a linear correlation between the price of the amnesiac lookback and its volatility, yielding a correlation coefficient is r <sup>2</sup> = 0.99. An interesting point is the difference in shape between the call and the put curves. As seen in **Figure 14**, the put curves slower than the call, due to the right vs. left tail of the lognormal distribution, especially when volatility is high in the case of ripple. This means that the amnesiac put may be useful if a trader anticipates a down market without needing to worry when to sell, while enjoying a greater range of values to choose from due to the slower convergence to the standard lookback option.

One can conceive of using Algorithm 1 to price an option, particularly for Ethereum which was developed for Smart Contracts via the Ethereum Virtual Machine (EVM). The scripting language it uses is very expressive and Turing complete, and the EVM acts as a distributed computer that executes commands based on the resource gas. Smart contracts thus become a good way to buy contracts that previously only existed in the OTC market. In this case, a transaction function transfers Ether to the writers account at the time of purchase. At the time of maturity, if the option is exercised, the smart contract executes the necessary transactions as determined by the option automatically. The buyer sends in a selection of monitoring periods A, then the input parameter could be the number of periods which is used to generate a subset of pseudorandom monitoring periods, or direct simulation of the selected periods. The latter produces a more accurate price but is also computationally more expensive and as a result, raises

TABLE 4 | Amnesiac call and put prices using historical volatility August 31, 2017 to November 30, 2017.


TABLE 5 | Amnesiac call and put prices using historical volatility December 01, 2017 to February 01, 2018.


the transaction fees of the smart contract. The choice becomes a question of cost vs. profit, and once statistical distributions on people's preferences of monitoring periods are collected, the pricing can be made even more accurate at lower computational costs.

## 5. CONCLUSION

Lookback options are immensely useful instruments for the use of hedging against risks associated with high volatility and notably effective at canceling investor regret as well. Although the literature has explored lookback options and many concepts of the partial lookback, previous work has mainly involved lookbacks with continuous monitoring, and equidistant lookback period selection. We suggest the use and exploration of a new form of partial lookback, the amnesiac lookback option, that can look back upon any discretely selected period prior to expiration.

We discuss the practical applications that amnesiac lookback options may have in the trading of cryptocurrencies such as Bitcoin and especially Ethereum. The cryptocurrency market has been defined by its high volatility trading and plethora of exogenous shocks that regularly disrupt the market. In such a high risk and speculative market, it is imperative that investors have some way of mitigating their risks and exposures. In this realm, lookback options are a simple and effective way to do so, while amnesiac lookback options allow an investor to fully customize the amount of risk he is willing to take on or forego by adjusting the number of lookback periods, and potentially where to place them, as shown in **Figure 1**. This is particularly true for put options under high volatility. Past studies have also found long memory in Bitcoin returns and return volatility, which path-dependent options such as the amnesiac lookback are able to correct by closing arbitrage opportunities, if traded on rationally. This opportunity is greatly enhanced by Smart Contracts, which decentralizes and automates highly customizable derivatives like path-dependent options.

Using Monte Carlo simulations under Black-Scholes assumptions, we find the Hill Function to be a suitable model for pricing the payoffs of the amnesiac lookback option. We also determine the bounds of the amnesiac lookback consistent with prior results on partial lookbacks, that is the standard lookback from above, and the corresponding vanilla European option (call or put) from below.

Our work additionally discusses the merits of equal spacing vs. random sampling of lookback periods. In this regard, we find that equally spaced lookback periods yield the greatest payoffs under Black-Scholes assumptions. Since the formula for modified trees and lattice methods produce analytic results for equidistant sampling, and equidistant sampling is a strict upper bound for random sampling, we can conclude that if days are not chosen equally then the instrument will be prone to being overpriced. More importantly, our results are more realistic and related to the designated function of these options in actual markets than some previous works.

The fact that how periods are selected affect the pricing suggest multiple paths of further inquiry. Since this is phenomenon under Black-Scholes assumptions, different models such as introducing stochastic volatility and Lévy processes may produce different results. Tying this to real commodities, statistical analysis of Cryptocurrency data can benchmark the efficiency of such an option. Simulating extreme, volatile events within the context of the Gaussian random walk would bridge the

empirical and axiomatic formulations, and is ultimately the purpose of having a choice in monitoring periods. Understanding period selection strategies will yield new perspectives of using path-dependent options for hedging risk, and expand trading strategies with the growth of novel and volatile blockchain commodities.

## AUTHOR CONTRIBUTIONS

H-CC wrote the code for the numerical simulation, data processing, and the mathematical finance literature review. KL contributed economic analyses, interpretation, literature review and extension to cryptocurrencies. Both contributed to the conception and design of the option's amnesiac feature, interpretation of results, and general authorship the paper.

## REFERENCES


#### FUNDING

Funding for this research was provided by the Jack Byrne Scholar Fund at the Department of Mathematics, Dartmouth College. Open access article fees were provided by the Dartmouth Open-Access Publication Equity Fund.

## ACKNOWLEDGMENTS

We'd like to thank Seema Nanda for proof-reading and overall suggestions, Feng Fu, Sergi Elizalde and Peter Winkler at the Department of Mathematics, Dartmouth College for their mathematical suggestions, Jennifer Kuo for statistical consulting and Jonathan Meng for web scraping consultation.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Chang and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# On the Aerodynamic Forces on a Baseball, With Applications

Gerardo J. Escalera Santos <sup>1</sup> \*, Mario A. Aguirre-López <sup>2</sup> , Orlando Díaz-Hernández <sup>1</sup> , Filiberto Hueyotl-Zahuantitla1,3, Javier Morales-Castillo<sup>4</sup> and F-Javier Almaguer <sup>2</sup>

<sup>1</sup> Facultad de Ciencias en Física y Matemáticas, Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Mexico, <sup>2</sup> Facultad de Ciencias Físico-Matemáticas, Universidad Autónoma de Nuevo León, San Nicolás de los Garza, Mexico, <sup>3</sup> Cátedra CONACyT, Mexico City, Mexico, <sup>4</sup> Facultad de Ingeniería Mecánica y Eléctrica, Universidad Autónoma de Nuevo León, San Nicolás de los Garza, Mexico

The aerodynamic forces acting on a baseball are those produced by the contact between the ball and the air, and are defined by the initial conditions of the pitch. It is well known that such forces determine the changes from the typical parabolic ballistic trajectory, either in the direction of the movement of the ball (drag force), or producing a lift or lateral deflection (Magnus and seam forces). The drag and Magnus effects have been widely studied and there are many references about their nature and the trajectory they produce, which is predictable. This has led to most baseball research being related with spinning pitches. On the other hand, there is the unpredictable motion of a knuckleball, whose erratic trajectory accompanied by a poor understanding of the forces produced by the asymmetry of the seams had markedly limited research about it until the beginning of this century. However, nowadays interest in the knuckleball is resurfacing. Data collected by wind tunnel experiments and real pitches have motivated researchers to analyze the phenomenon and build models that try to predict the motion of the ball. In this work we aim to provide the reader some basic ideas on aerodynamic forces through a combination of experimental results, phenomenological and dimensional analysis, with special focus on new advances on the seam effects of a knuckleball pitch. In addition, we discuss possible ways to extend the existing models about the seam forces. Finally, we summarize from the literature some methods regarding the reproduction and reconstruction of baseball trajectories from aerodynamic forces and discuss their application as well.

Keywords: baseball, knuckleball, seams, drag force, lift force, Magnus force, ball games

## 1. INTRODUCTION

Aerodynamics and the flight of baseballs are very interesting phenomena from the point of view of sports, engineering and science. The different types of pitches are denoted by the different initial configurations that the pitcher gives to the ball by means of his hand. Each initial combination of velocity and spin produces a way of interaction between the air and the ball which results in curveballs, fastballs, sliders, and knuckleballs, among others. Moreover, such interaction can be affected by other factors, such as when the ball is dented or when it has liquid on its surface. This can cause an erratic movement in the expected trajectory of the ball, and be the difference between a strike or a home run.

#### Edited by:

Jian-Ao Lian, Texas A&M University System, United States

#### Reviewed by:

Kazuharu Bamba, Fukushima University, Japan Yonghui Wang, Prairie View A&M University, United States

#### \*Correspondence:

Gerardo J. Escalera Santos gescalera.santos@gmail.com

#### Specialty section:

This article was submitted to Mathematics of Computation and Data Science, a section of the journal Frontiers in Applied Mathematics and Statistics

> Received: 28 March 2018 Accepted: 21 December 2018 Published: 28 January 2019

#### Citation:

Escalera Santos GJ, Aguirre-López MA, Díaz-Hernández O, Hueyotl-Zahuantitla F, Morales-Castillo J and Almaguer F-J (2019) On the Aerodynamic Forces on a Baseball, With Applications. Front. Appl. Math. Stat. 4:66. doi: 10.3389/fams.2018.00066

For engineering, the evolution of a baseball during its flight means an open door to many possible applications. In turn, for science it means the study of the process in which a subsonic flow interacts with a solid rough sphere, having the characteristic that for zero or lower spin values of the spin the trajectory is very erratic and becomes unpredictable, whereas for higher values the trajectories are smooth and predictable. In this way, the trajectories of a baseball are commonly classified by spinning pitches and non-spinning pitches.

The predictability of both types of trajectories seems contradictory at first glance, since one would expect the spinning balls to have a more complex movement because more forces act on them. This is not the case, however, and the main cause of this is the role of the seams of the baseball. When the interaction between the seams and the air is insignificant, more predictable and stable trajectories occur, and vice-versa. This last statement can be better understood by looking at the Reynolds' transport theorem, which establishes [1–3]

$$\frac{\partial}{\partial t} \int\_{\mathcal{V}} \rho \mathbf{V} d\mathcal{V} + \int\_{\mathcal{S}} \rho \mathbf{V} \left(\mathbf{V} \cdot \mathbf{n}\right) dA = \int\_{\mathcal{S}} \mathbf{T} \cdot \mathbf{n} d\mathcal{S} + \int\_{\mathcal{V}} \rho \mathbf{b} d\mathcal{V} \text{ (1)}$$

for a fluid-containing volume of space with volume V and surface S, where **V** and ρ denote the velocity and density of the fluid, respectively, and **n** is a unit normal vector pointing outside S. In turn, the right-hand side denotes the sum of all forces acting on the volume of space, which can be classified according to the intensive properties as body (first term) and surface (second term) forces, with **b** and **T** being the body forces vector per unit mass and the stress tensor, respectively. Among body force gravity, Coriolis and centrifugal forces can be found, whereas the forces generated at the ball-air boundary include pressure, normal and shear stress, etc.

In spite of gravity, body forces are weak and are commonly neglected in calculations of ball sports as reported in calculations by Robinson and Robinson [4] and Aguirre-López et al. [5]. On the other hand, the surface forces can be classified according to the direction and nature of their origin: drag, Magnus, lift, side and other forces [6–8]. Magnus force can be defined as that caused by the spin of the ball, therefore it is present only in spinning pitches. In turn, we define lift and side forces as those caused only by the motion of the baseball through the air without rotation (non-spinning pitches); then it depends on the orientation of the ball because of the asymmetry of the seams. A proposed curve of how lift and Magnus forces behave when varying the magnitude of the spin is shown in **Figure 1**. The lift and side forces produced by the seams decrease with the increase of the angular velocity (which is expected according to Watts and Sawyer [11], Mehta [12], and Cross [9]) while by definition, the Magnus force approaches zero when the angular velocity vanishes [8]. In this way, the non-spinning window is related to erratic and unstable trajectories whereas the spinning window is related to predictable and stable trajectories. Finally, there is a window between both cases, in which the superposition of forces is more significant and the ball can continue with the erratic movements of non-spinning pitches or can draw a smooth, slow

FIGURE 1 | Scheme of the maximum value of lift, side and Magnus forces when varying the angular velocity ω for a given speed V. The y−axis is scaled according to the maximum drag force for the considered velocity. The curves for lift, side and Magnus forces are drawn according to experimental data [8–10]. The missing curve for lift and side forces between spinning and non-spinning windows is because there are no data reported in such range.

and predictable trajectory. From our understanding, this is the main reason why pitchers do not throw balls within this range.

The curves in **Figure 1** also show the evolution of the study of the types of pitches at this time. Until the beginning of this century, spinning pitches covered most baseball research; therefore, books and compilations of works related to the subject avoided the phenomena present in non-spinning pitches [7]. Nowadays, the interest in non-spinning pitches is resurfacing and there is more information available in books of the present decade, like Cross [6]. However, the new developments and methods for studying non-spinning and intermediate windows are very disconnected, so a new compilation of work on the matter is necessary.

The purpose of this work is to introduce the reader to the diverse existing methodologies for studying the aerodynamic forces in the flight of a baseball, and to discuss the possible ways to complement such studies for future research and applications. All this is based on a literature review of the subject. Each topic begins with a brief derivation of the mathematical model of the respective force; then the model is compared with the results obtained from experimental measures, and a discussion of this is achieved. We have special interest in discussing how the seams affect the aerodynamics of the ball; therefore, the lift force and non-spinning pitches will get more attention.

This review is structured according to the classification of pitches mentioned above. In this way, we begin by presenting the most common mathematical model to describe the movement of a baseball with a concise review of the spinning pitches in section 2. Next, section 3 deals fully with the advances on nonspinning pitches and presents some comments on the boundary layer and the wake of a baseball in motion. Finally, we present a compilation of the existing uses regarding the aerodynamics of baseballs, and outline the research trends of the subject in section 4.

## 2. SPINNING PITCHES

In accordance with **Figure 1**, we define a spinning pitch as any throw, excepting knuckleballs, that has an initial velocity in the range **V** = [(Vx,Vy,Vz)] ∈ [(−3, 30, −3), (3, 50, 3)] m/s and an initial spin of |ω| = [100, 310] rad/s, which are the values at which baseballs are thrown by professional pitchers when fixing the y−axis in mound-home direction and z−axis perpendicular to the floor, according to the right-hand rule [13]. In this way, some examples of spinning pitches are the curveball, slider, change-up, and the variants of the fastball. The hallmark of this type of pitches is the Magnus force. We begin by describing the drag force in section 2.1, which is present in all types of pitches, with the aim of establishing some basic concepts to introduce the Magnus force.

#### 2.1. The Drag Force

Drag or friction is the force that resists the movement of an arbitrary object. The classical way to derive drag force is by considering that the study of a ball moving through a static medium is equivalent to the one of a static ball with fluid in motion; therefore, **V** corresponds to the velocity of the fluid, and thus a mathematical formula for the drag of the ball can be obtained from the momentum conservation equation [2, 14]

$$\frac{\partial \mathbf{V}}{\partial t} + (\mathbf{V} \cdot \nabla)\mathbf{V} = -\nabla \left(\frac{\rho}{\rho}\right) + \mathbf{g} \qquad , \rho \equiv \text{pressure}, \quad \text{(2)}$$

which, in turn, is computed by considering that the air is a Newtonian and incompressible fluid in the Reynold's transport theorem (1); see Ferziger 1996 [2] for the computing of Equation (2). The conservation of mass equation also reduces to

$$\nabla \cdot \mathbf{V} = 0.\tag{3}$$

Then, assuming a steady flow the Bernoulli's equation results

$$\frac{1}{2}V^2 + \frac{p}{\rho} \equiv \text{constant},\tag{4}$$

where we used <sup>1</sup> 2 ∇V <sup>2</sup> = **V** × (∇ × **V**) + (**V** · ∇)**V** avoiding the gravitational term [14]. From Bernoulli's Equation (4) the reader can observe that the greatest pressure occurs at points where the velocity equals zero. These points are commonly found on the surface of solid bodies and are usually called stagnation points. For the case of a sphere ball, one stagnation point is present on the front side. The pressure at this point is

$$p\_{\text{max}} = p\_0 + \frac{1}{2}\rho V^2 \tag{5}$$

where p<sup>0</sup> is the pressure of the fluid at infinity [14, 15]. In this way, the second term is related to the force opposing the motion of the ball, the drag force, and thus drag can be calculated by the difference of pressure between the front and back sides of the ball [2]. However, there are other phenomena affecting the drag. At the corresponding Reynolds number (Re) of a typical throw (10<sup>4</sup> − 10<sup>5</sup> ) the flow is not steady but turbulent [8, 10, 13, 16, 17]; then the point of separation of the boundary layer moves away from the front of the ball when the ball's velocity increases, making the wake smaller and avoiding a momentum in the reverse direction on the rear of the ball [18, 19].

The approximation of the drag force is done by introducing the sectional transverse area A and a dimensionless coefficient C<sup>d</sup> into the second term of the right-hand side in the Equation (5), such that

$$\mathbf{F}\_d = -\frac{1}{2}\rho \mathbf{C}\_d AV \mathbf{V}.\tag{6}$$

where C<sup>d</sup> = Cd(V), and then C<sup>d</sup> acts in (6) as a fraction of the area A interacting with the air, or in other words, a measure of how aerodynamic the ball is [18, 20, 21].

As a consequence, the study of the drag coefficient (Cd) is so important that it becomes the most effective formula to approximate the drag; therefore, the drag force (**F**d) is commonly expressed in terms of Cd. There are two ways to explain the behavior of Cd. On the one hand, there is the drag crisis phenomenon for smooth spheres mentioned in Landau and Lifshitz [14], which indicates that a crisis in the value of C<sup>d</sup> can occur at Re∼ 10<sup>5</sup> (at velocities of 30–40 m/s for a typical baseball). Data shown in Frohlich [22] and Cross [23] fit to such model. On the other hand, from the famous curve of Adair [7] to the recent data obtained by Naito [24], most of the experimental measures suggest a model without drag crisis [6, 17, 25, 26]. This predominance maintains for other sphere balls [23, 27–29] and gives rise to most of the used models to compute the drag force without including the drag crisis phenomenon. There are two models of special interest: the model by Cross [6], from which the drag coefficient can be obtained if the instantaneous ball's velocity at a fixed time is known, and the curve of Adair [7], which approaches to

$$C\_d(V) = 0.29 + 0.22 \left[ 1 + e^{(V - 32.37)/5.2} \right]^{-1},\tag{7}$$

according to Aguirre-López etal. [5]; see **Figure 2**.

Finally, it is important to mention that the drag has an oscillating dependence on the orientation of the seams. However, for spinning pitches such oscillations average to zero and so they are not considered here. The dependence of C<sup>d</sup> on the seams' orientation will be discussed in section 3.1.

#### 2.2. The Magnus Force

Magnus force is the essential characteristic of a spinning pitch, as we mentioned before. The Magnus effect is observed as smooth deflections in the trajectory of a ball. All of us have a clear empirical knowledge of the Magnus effect: large deflections are reached by increasing the spin frequency of the ball. Although this statement is true, the direction of the deflection varies for different configurations of linear velocity **V**, the angular velocity ω and the ball properties. In fact, for any viewer, the expected direction of the ball's deflection goes on (ω × **V**),

as illustrated in **Figure 3**. However, a reverse direction of the Magnus force has been reported for smooth balls like those used in soccer [27], and also in smooth spheres simulating baseballs in experiments by Briggs [30], Cross and Lindsey [31]. This phenomenon is commonly called the "anti-Magnus effect" and it is possible only for a range of Re and spin when one side of the smooth ball remains in a laminar flow while the opposite side becomes turbulent. Then, a low pressure region is originated in the turbulent side because it is generally farther to the ball surface than the laminar layer. Thus, the ball moves to the region with lower pressure by conservation of momentum, as illustrated in **Figure 3**. For a detailed explanation of the causes of the reverse in the direction of Magnus force, the reader is referred to Cross and Lindsey [31]. For a general understanding of the phenomenon, the reverse-Magnus occurs at Re∼ 10<sup>5</sup> combined with a low speed due to Magnus force Rω (where R is the radius of the ball) compared with the ball's velocity V, such that the spin factor S = Rω/V is in the range 0–0.6.

However, the baseball is not smooth and there are no reported studies of an anti-Magnus effect in baseballs. This is because the seams of the ball give to it some roughness that accelerates the separation of the boundary layer in both the up and down sides of the ball. Indeed, the Magnus effect in baseballs' flight arises because one side of the ball offers larger friction than its opposite side, which means that the speed of the main flow of air is larger on the former and as a consequence the lower-pressure region locates at the opposite side, according to the anti-Magnus phenomenon [6–8]; see the schemes in **Figure 3**. Therefore, the nature of the Magnus force is similar to that of the drag force, due to a difference of pressure. For this reason, Magnus force is commonly written in a similar way as drag force (Equation 6), namely,

$$F\_{\rm M} = \frac{1}{2}\rho A C\_{\rm M} V^2,\tag{8}$$

for an arbitrary direction of motion, with C<sup>M</sup> being the Magnus coefficient. Moreover, in the same system of coordinates to drag force, Equation (8) becomes

$$\mathbf{F\_M} = \frac{1}{2}\rho A C\_\mathbf{M} \sin\phi \, V^2 \hat{\mathbf{n}},\tag{9}$$

where the unit vector **n**ˆ = ω×**V** |ω×**V**| gives the direction of the resulting linear momentum, φ is the angle between ω and **V** so that (V sin φ) is the component of V that contributes to the force [6, 8, 10]. It is important to remark that equation (9) has been widely used to reproduce Magnus force of sport balls and other areas of aerodynamics [32–34], and it has been proved experimentally (for sport balls) only for φ = 0, 90 and 180◦ [8, 10].

The Magnus coefficient C<sup>M</sup> is a function of ω and V for arbitrary magnitudes of such variables because Equation (9) depends on both the instantaneous and spin velocities [6, 10, 31]. However, Nathan [17] has reported that C<sup>M</sup> behaves independently of V for ω and **V** values in the range at which spinning pitches are thrown. Moreover, the value of C<sup>M</sup> in such range is similar to the Magnus coefficient for other sport balls, and there are some models to calculate its value [8, 30, 35]. **Figure 3** shows the model (10):

$$\mathcal{C}\_{\rm M}(\omega) = 3.19 \times 10^{-1} \left[ 1 - e^{-2.48 \times 10^{-3} \omega} \right],\tag{10}$$

which has been used to compute the Magnus force in numerical simulations [5, 8].

#### 2.3. Discussion and Potential Research Trends

In the following items, we summarize and discuss the highlights of the spinning pitches, and also sketch out the potential research trends of the subject:


Both of them could be promising areas of opportunity to characterize the drag, and the aerodynamics of a baseball, in a more comprehensive way.


## 3. NON-SPINNING PITCHES

Non-spinning pitches consist of a specific type of throw: the knuckleball. It contains all combinations for the linear and angular velocities in the ranges V ∈ [20, 40] m/s and ω ∈ [0, 50] rad/s, respectively [36]. Knuckleball pitches are the most interesting throws for aerodynamics because the ball can no longer be considered a sphere since the effect of the seams is significant. However, this is more complex to understand and, therefore, knuckleball studies are fewer than those for spinning pitches.

In this section we discuss the advances in drag force for non-spinning pitches (section 3.1). Then, the description of the models for lift and side forces is presented in section 3.2. An introduction to the modeling of the seams and the boundary layer observations is shown in section 3.3. A compilation of the knuckleball's research and trends is discussed in section 3.4.

## 3.1. The Drag Force in Non-spinning Pitches

Drag force is different in knuckleballs than in spinning pitches. The drag is approximately constant for a specific velocity in a spinning ball, however, in a knuckleball it is not. As mentioned in section 2.1, the Re of a baseball pitch corresponds to an unsteady flow and then an experimental coefficient must be introduced in the model for the drag force (6). However, when the ball does not spin some vortices are shed from the ball and then the drag oscillates in time. Ferziger and Peric [ ´ 2] discusses such an effect for a smooth cylinder simulated by CFD. The drag on the cylinder oscillates periodically with a frequency according to the appearance of the vortices such that it has one maximum and one minimum during the formation and shedding of each vortex. Such vortex shedding has also been reported for baseballs in Texier et al. [37], whose effect in lift and side forces will be discussed in section 3.2.

In addition to the variation in time, the drag changes when varying the orientation of the ball. Investigations carried out in the present decade show structured oscillations of drag coefficient despite the turbulence present in the phenomenon. For example, the experiment by Higuchi and Kiura [38] with different configurations of the ball, namely, the four-seam (4S), the twoseam (2S) and an arbitrary orientation of the ball<sup>1</sup> . They found the largest variation in oscillations for the 4S orientation, which is about twice as large as the case of the 2S orientation and around four times that of the arbitrary orientation. The shape of the oscillations in the 4S orientation is maintained for ball velocities in the range of 16–30 m/s.

Similar average variations for C<sup>d</sup> have been reported by Alam et al. [25] in studies of the drag force for Major League Baseball (MLB) and National Collegiate Athletic Association (NCAA) baseballs and softballs. They reported lower variations in C<sup>d</sup> values for NCAA than for MLB baseballs at different ball orientations, which suggests a dependence on the height of the seams (1.5 mm for NCAA and 1mm for MLB balls). At first

<sup>1</sup>A detailed explanation of the typical baseball orientations can be found in Borg and Morrisey [36].

glance, the idea seems to be solid because a larger "extra-obstacle" should cause a corresponding "extra-drag" in the fluid. Even more, this is supported by the results of Kensrud et al. [39], who analyzes hit balls with different heights of seams and found that balls with smaller seams reach larger distances, which indicates a lower drag.

Finally, it is important to remark that commonly C<sup>d</sup> decreases with increasing V for baseballs, softballs, cricket balls and smooth spheres. In the case of baseballs, the value of C<sup>d</sup> decreases from ∼ 0.6 to ∼ 0.4 dimensionless units [25, 32].

#### 3.2. The Lift and Side Forces

The main reason for which knuckleball pitches are much more unpredictable than spinning pitches is because of lift and side forces<sup>2</sup> . Similar to the drag and Magnus forces, lift force is produced by an imbalance in pressure and then it is proportional in magnitude to the square of the ball velocity [11, 40], so that

$$F\_L = \frac{1}{2}\rho C\_L A V^2. \tag{11}$$

However, the behavior of the lift coefficient (CL) is not like C<sup>d</sup> or C<sup>M</sup> but it varies with the angle of attack θ of the ball. **Figure 4** shows how unpredictable a knuckleball can be, even for the most symmetric ball orientations (4S and 2S). Results of Borg and Morrisey [36] show four cycles in C<sup>L</sup> for the 4S orientation, each one with a period of 90◦ , with a semi-sinusoidal behavior and ringlets at the end of a cycle. This means that a variation of ∼ 22.5◦ in smooth-angle zones may or may not produce a maximum/minimum lift; instead, a variation of only ∼ 10◦ in the zone of ringlets can produce any type of motion. In addition to such complexity for a strictly non-spinning ball, in practice it is difficult for a ball to travel without rotating since a little spin is induced by contact from the air with the seams [12, 37]. As a consequence, the ball's trajectory could have an apparent erratic motion whereas the map of balls passing through the home plate in a real pitch is seen as random when varying θ [40]. All of this makes it difficult to compute a model that may accurately reproduce the trajectory of a knuckleball with any orientation.

It is important to mention that there is controversy in the causes that originate the lift force. It is evident that the asymmetry of the seams plays a fundamental role in causing such force. According to Watts and Sawyer [11] and Mehta [12], there are two possible ways to produce a lift force on a baseball: for many years, the classical hypothesis stated that the lift is produced not only by the seams but by the shedding of vortices that occurs at the rear side of the ball. All of these interact in a complex way, as mentioned in Ferziger and Peric [ ´ 2], Watts and Sawyer [11], and Mehta [12]. On the other hand, and according to Texier et al. [37], there is the possibility that the lift force may be originated only by the perturbations at the front side of the ball. This means that the seams produce the total lift of the ball, especially those located at the separation or critical points of the boundary layer at 52, 140, 220, and 310◦ [12, 36]; this will be addressed in detail in section 3.3. Texier calculated that the force produced by the vortices at Strohual numbers (St) of St∼ 0.2 is significantly lower than the magnitude of the lift measured in experiments, so that it practically does not contribute to lift. More information about lift force in ball sports can be found in investigations by Mehta [12], Hong et al. [42], and Murakami et al. [43].

#### 3.3. The Seams and the Boundary Layer

The manner in which the seams affect the boundary layer is very sensitive to small changes between seams. As commented by Borg and Morrisey [36], when a seam is located near the natural separation angle (the angle of separation of a smooth ball), it can induce turbulence and consequently provoke a delay in the separation of the boundary layer; in turn, the seam can force the separation to occur and cause an advance in the separation angle. Such effects are seen in the sudden changes in the values of the lift coefficient for a 4S orientation. For a smooth ball, a sinusoidal shape of C<sup>L</sup> would be expected, in a similar way to the lift coefficient reported for a smooth cylinder by Ferziger and Peric [ ´ 2]. However, **Figure 4** suggests a quasi-periodic behavior for C<sup>L</sup> for 4S balls, with fast changes in magnitude and direction at 52, 140, 220 and 310◦ , as observed by Watts and Sawyer [11]. This is because the separation point is located around such degrees and then it advances or delays with a little variation in the stitch position.

Experiments by Higuchi and Kiura [38] show that a variation of only one degree (36 to 37◦ ) in the stitch position causes a sudden separation. Moreover, they reported that the balls are more susceptible to hysteresis (including induced rotation) at the zones of separation. For 4S balls and Re above 1.5×10<sup>5</sup> , the ball is sensitive to the initial rotation, namely, spins of 5 rad/s become 10.5 rad/s, increasing linearly and having a spin limit of 18.9 rad/s even for Re above 2×10<sup>5</sup> . In turn, for 2S balls they found that the oscillation frequency is constant over Re∈ [1.9 × 10<sup>5</sup> , 4.6 × 10<sup>5</sup> ]. As a consequence of the induced rotation, the phenomenon becomes more unpredictable because the separation point moves forward or backward at every moment of time. This is the reason why the throws inside the intermediate window in **Figure 1** are the most difficult to study. We invite the reader to consult the research of Higuchi and Kiura [38] for detailed observations of the boundary layer.

To end the collection of the advances on knuckleballs, it is important to mention the phenomenological model proposed by Aguirre-López et al. [41] for computing the lift coefficient. It consists of computing a super-imposition of the forces produced by the vortex shedding and each stitch, so that

$$\mathcal{C}\_{L}(\theta) = a\_{0} \sin(4\theta - \pi) + a\_{1} \sum\_{i=1}^{n} \left[ \sin \left( \frac{||\mathbf{s}\_{i} - \mathbf{p}||\pi}{2d} + \pi/2 \right) \right] \tag{12}$$

$$\text{sgn} \left( \rho^{\*} - s\_{i}^{\*} \right) \Big|, \tag{12}$$

where C<sup>L</sup> = CL(θ) is now a function of the angle of attack of the ball, the first term in the right-hand side is the force caused by the vortex shedding and the term P(·) is the sum of forces produced by the seams, **p** is the stagnation point, **s**<sup>i</sup> is the position of the i-th stitch, s ∗ i and p ∗ are the z−components of **s**<sup>i</sup> and **p**, respectively, and a<sup>0</sup> and a<sup>1</sup> are weight coefficients.

<sup>2</sup>From this point on we will use the term "lift force" for referring to both lift and side forces because both are equivalent, excluding gravity.

Therefore, each stitch in model (12) produces a force whose magnitude decreases smoothly when the stitch moves away from the stagnation point and takes into account the symmetry on the z-axis by introducing the sign function sgn(·). Thus, despite the fact that model (12) does not consider effects of hysteresis and high sensitivity to perturbations at the separation point, it fits the experimental data of Borg and Morrisey [36] for 4S and 2S orientations, as shown in **Figure 4**. Model (12) opens the door to future research on how the seams and the vortex shedding affect the lift force.

## 3.4. Discussion and Potential Research Trends

Here we summarize and discuss the highlights of the nonspinning pitches and outline potential research trends of knuckleballs pitches as follows:


## 4. APPLICATIONS

As the reader may suppose, there are many ways to make useful the information compiled in sections 2 and 3. And, indeed, baseball studies have been the basis of numerous technologies on the matter, specifically those ones about spinning pitches. We finalize this work with a brief summary of the main applications of the aerodynamic forces on baseballs. Section 4.1 is focused on the studies related to baseball's trajectories. In turn, section 4.2 talks about the complementary applications including the best known of them: the PITCHf/x algorithm.

## 4.1. Prediction and Reconstruction of Trajectories

#### 4.1.1. The Simulation Problem

The most simple use of aerodynamic forces is the simulation or prediction of trajectories. A simulation of a baseball pitch is frequently carried out by using the Equations (6) and (9) along with gravity to compute a model of forces as:

$$m\dot{\mathbf{V}} = \mathbf{F}\_d + \mathbf{F}\_M + m\mathbf{g},\tag{13}$$

which can be solved numerically by Runge-Kutta-4 or other integration methods, whereas drag and Magnus coefficients can be computed by Equations (7) and (10) or similar approximations [6, 44]. The simulations could become more realistic by including eventual forces in the model (13). An example of this is the model of Robinson and Robinson [8], which adds a constant in wind to the ball velocity so that **V** ′ is redefined as (**V** ′ = **V** + **W**), where **W** is the wind velocity.

In turn, simulations of baseball trajectories are commonly applied to some sport and technology areas such as in video games [45], baseball machines [10] as well as for instruction for baseball players. The last one is the main reason for which research on knuckleballs is a topic of special interest.

#### 4.1.2. The Reconstruction Problem

The counterpart of simulation is the reconstruction of trajectories. In these works, there is no possibility to give the initial conditions of a pitch and obtain the trajectory but the trajectory must be extracted, tracking or reconstructed only by a set of 3D or 4D (space plus time) points belonging to the trajectory. Such data points are commonly recorded by baseball broadcast videos and/or images of real games; therefore, the trajectory obtained is used for the replay in television.

Various methodologies have been reported for extracting or tracking the trajectory. Most of them use types of diverse filters like color, position, size and shape [46–48], and others [49, 50] select possible trajectories. Then the chosen trajectories are compared with the model (13) so that, if the resulting trajectory does not agree with the model, then it is discarded and a new one is needed. Takahashi et al.'s investigation [47] also deals with classifying the type of pitch by relating a total of 36 features, including the shape and speed calculated from the ball trajectory data and the ball speed from the screen display. They report an accuracy of ∼ 89% with their methodology.

The methodologies that deal with a "direct" reconstruction are based on the use of the equations of motion. Shum and Komura [51] and Miyata et al. [52] use color filtering for detecting 2D candidate trajectories. Then, Shum and Komura estimate the depth of the ball in the scene by introducing a model (13). In turn, Miyata et al. [52] chose one candidate by fitting a uniformly accelerated motion model [similar to model (13)] and finally, they use multiple cameras calibrated temporally and geometrically to obtain a 3D trajectory. On the other hand, Aguirre-López et al. [5] developed an algorithm that directly solves the model (13) in two interrelated parts by decoupling the Magnus force from the equations of motion, using the Newton-Raphson method when knowing three points of the trajectory. They reported absolute error values of ∼ 0.1 mm between simulated and reconstructed trajectories. Finally, Kagan and Nathan [53] have developed a software called the trajectory calculator, which is similar in operation to that of Aguirre-López et al. [5] but with simpler assumptions. The results are less accurate but it is a good tool to start in the subject. The trajectory calculator can be downloaded directly at [54].

#### REFERENCES


#### 4.2. The PITCHf/x Algorithm and Clustering

The second part of the application deals with problems related to the classification of trajectories, among which the PITCHf/x algorithm is the most popular and accurate reported method in research and in the world of baseball. The algorithm (including new versions and software packages) has been consolidated as a powerful tool in the area of pitch classifications [55, 56].

The PITCHf/x algorithm consists of two parts. The first one involves reconstructing trajectories by estimating the coefficients Cd, C<sup>M</sup> and the spin axis φ (the angle between y−axis and ω) using non-linear least-squares fitting with the Levenberg-Marquardt algorithm. Nathan [55] reports very good adjustments; indeed, root-mean-square deviations of the fitted trajectory of around 1 mm in each dimension. The second (and the main) part of the work deals with the classification of pitches. The classification is based on (**V** vs. φ) and (ω vs. φ) graphics. As a result, the types of pitches are arranged in clusters in polar scatter plots and scatter plots of the deflection of the ball at home. Pane [57] carried out an interesting cluster analysis from the results of Nathan [55] based on PITCHf/x. The research on the topic continues.

#### ETHICS STATEMENT

We express our gratitude to Elsevier for the permission to reuse a modified version of the copyrighted **Figures 2**, **3a**, **4** in the present paper. Order numbers: 4511461494464 and 4511510076028.

## AUTHOR CONTRIBUTIONS

GE and OD-H contributed to the design and discussion of this document. MA-L contributed to the design and writing of this document. FH-Z contributed to the discussion and writing of this document. JM-C and F-JA contributed to the discussion of this document and as advisor.

#### ACKNOWLEDGMENTS

FH-Z thanks CONACyT, Cátedra 873.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Escalera Santos, Aguirre-López, Díaz-Hernández, Hueyotl-Zahuantitla, Morales-Castillo and Almaguer. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.