Great Barrier Reef degradation, sea surface temperatures, and atmospheric CO2 levels collectively exhibit a stochastic process with memory

Quantifying ecological memory could be done at several levels from the rate of physiological changes in an ecosystem all the way down to responses at the genetic level. One way of unlocking the information encoded in a collective environmental memory is to examine the recorded time-series data generated by different components of an ecosystem. In this paper, we probe into the case of the Great Barrier Reef (GBR) which is threatened by elevated sea surface temperatures (SST) and ocean acidification attributed to rising atmospheric CO2 levels. Specifically, we investigate the interrelated dynamics between the degradation of the GBR, SST, and rising atmospheric CO2 levels, by considering three datasets: (a) the mean percentage hard coral cover of the GBR from the archives of the Australian Institute of Marine Science; (b) SST close to the GBR from the National Oceanic and Atmospheric Administration; and (c) the Keeling curve for atmospheric CO2 levels measured by the Mauna Loa Observatory. We show that fluctuating observables in these datasets have the same memory behavior described by a non-Markovian stochastic process. All three datasets show a good match between empirical and analytical mean square deviation. An explicit analytical form for the corresponding probability density function is obtained which obeys a modified diffusion equation with a time dependent diffusion coefficient. This study provides a new perspective on the similarities of and interaction between the GBR’s declining hard coral cover, SST, and rising atmospheric CO2 levels by putting all three systems into one unified framework indexed by a memory parameter μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} and a characteristic frequency ν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}. The short-time dynamics of CO2 levels and SST fall in the superdiffusive regime, while the GBR exhibits hyperballistic fluctuation in percent coral cover with the highest values for μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} and ν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}.


Introduction
Ecological memory where past events shape and alter future evolution of ecosystems in a changing environment has been a recurring theme in understanding the response of an ecological community to climate change (Hughes 2019;Schweiger et al. 2018;Peterson 2002;Ogle 2015;Nyström and Folke 2001;Johnstone 2016). An investigation of the Great Barrier Reef (GBR) demonstrated such memory in the mass coral bleaching that occurred in 2016 and 2017 (Hughes et al. 2019). Memory can play an ancillary role as an ecosystem acclimatizes and develops resilience when subjected to recurring climate extremes where adaptation for coral reefs can occur even at the genetic level (Thomas and Palumbi 2017;Romero-Torres et al. 2020;Carballo-Bolaños et al. 2020;Mumby and van Woesik 2014;Sully et al. 2019;Jurriaans and Hoogenboom 2019;Fine et al. 2019;Davidson 2019). The importance of climate memory was, in fact, underscored as being crucial in determining tipping points or critical transitions in climate-sensitive systems (van der Bolt et al. 2018;Turner et al. 2020). Oceans, as grand matrix of ecosystems, have a perceived slow dynamics and long memory (Payne et al. 2017). A manifestation of environmental memory of a reef coral was also observed in a 10-year study in Phuket, Thailand (Brown et al. 2015). Much earlier, ecological memory was investigated in relation to a reef matrix functional diversity to maintain itself following anthropogenic disturbances (Nyström and Folke 2001). Climate predictions also depend on climate memory embodied in multi-year recorded events that capture the past needed in fine-tuning computational models to improve forecasting reliability (Corti et al. 2012;Davini et al. 2017;Yuan et al. 2019). In this paper, we use a stochastic process with memory to provide a common analytical framework in capturing ecological memory as an interconnected phenomenon between the rising atmospheric CO 2 levels, an elevated sea surface temperature (SST), and the observed GBR degradation. Decades of continuing rise in atmospheric CO 2 with direct impact on rising sea surface temperatures that coral reefs must continually adapt to for survival give rise to a collective cause and effect memory. Putting atmospheric CO 2 , sea surface temperatures, and GBR degradation in one unified stochastic framework gives us a better comparative scale and gauge of this interconnected collective ecological memory.
Components of an ecosystem may operate at different time scales, but they could reveal a behavioral imprint embedded in a collective memory of the rising atmospheric CO 2 that elevates sea surface temperatures leading to the GBR degradation. Extensive studies have shown a consistent heat gain by the earth system for the past decades. Energy, in the form of heat, leaving and entering the earth climate system is continually measured in terms of radiation at the top of the atmosphere (von Schuckmann 2020; Meyssignac et al. 2019). Solar radiation heats the earth which in turn radiates back to space the reflected short wavelength radiation and longer wavelength infrared radiation (von Schuckmann 2020; Liu et al. 2020;Meyssignac et al. 2019). Re-emitted infrared from earth, however, is absorbed by atmospheric CO 2 and water vapor that cover the planet giving rise to a natural greenhouse effect. Hence, sunlight absorbed by earth is on average greater than outgoing longwave radiation (Carballo-Bolaños et al. 2020;Meyssignac 2019;Anderson et al. 2016). Greenhouse gases such as CO 2 , which trap infrared radiation, have been continually increasing in the atmosphere due to fossil fuel combustion and other human activities (Hughes et al. 2017a, b;Hansen 2013;Yue and Gao 2018;Canadell et al. 2007;Corti et al. 1999). The energy accumulating in the climate system for decades has inevitably contributed to raised sea surface temperatures.
Using satellite and ocean data, a thorough study on how energy is accumulated and distributed within the climate system for the period 1985-2017 was recently done . The biggest thermal reservoir on earth are the global oceans that cover most of the earth's surface with high heat capacity and absorbing more than 90% of earth's accumulated energy (Meyssignac 2019;Abraham et al. 2013;Cheng et al. 2017;Cheng et al. 2019;von Schuckmann 2020). Since 1998, examined ocean basins have exhibited increased warming not only in the top ocean layer but also at deeper layers of the ocean with increasing ocean heat content (Cheng et al. 2017;Liu et al. 2020;Frade 2018;Schramek et al. 2018) whose long-term effects could impact marine ecosystems such as the GBR, which stretches more than 2300 km, being visible from outer space. The GBR has been losing more than 50% of its initial coral cover for the past decades (De'ath et al. 2012;Mongin 2016;Hughes et al. 2018;Hughes et al. 2019).
Rising levels of atmospheric CO 2 , in fact, poses a double threat to the GBR. Aside from elevated sea surface temperatures, ocean acidification also occurs when excess CO 2 in the atmosphere is absorbed by seawater. When CO 2 enters and dissolves in seawater, it produces a chemical reaction that yields H + ion concentration effectively lowering the pH of the ocean (Brown et al. 2019;Albright et al. 2016;De'ath et al. 2012;Hughes et al. 2017a, b;Mongin et al. 2016;Osborne 2017;Orr 2005;Anthony et al. 2011). The GBR and levels of atmospheric CO 2 , together with SST, constitute complex systems subjected to a number of intractable factors. The observables corresponding to these systems are characterized by fluctuations arising from short-term and long-term influences in time and in geographical location which may be periodic or non-periodic (Myers et al. 2018) exhibiting a collective ecological memory. The concern over the fate of the GBR has spurred various ways of modeling coral reef dynamics ranging from hydrodynamic models (Lambrechts et al. 2008;Mongin et al. 2016) to solving differential equations (Melbourne- Thomas et al. 2010;Wolanski et al. 2004), numerical models (Andutta et al. 2013;Zychaluk et al. 2012), and stochastic processes with varying assumptions (Andutta et al. 2013). Nonlinear correlations with past states may exist rendering Markovian modeling for time series data sets insufficient. Modeling fluctuating observables in the interrelated dynamics of the coral cover of the GBR, the rise in atmospheric CO 2 levels, and SST, therefore, finds a non-Markovian stochastic process most appropriate.
In climate models, introducing randomness or stochasticity is crucial in representing the many intractable unresolved factors that determine the fate of interconnected ecological systems. Stochasticity has been incorporated, for example, in the project Climate SPHINX (Stochastic Physics High Resolution Experiments) to predict future climate scenarios (Davini et al. 2017;Meccia et al. 2020;Strommen et al. 2019;Yang et al. 2019). The inclusion of stochastic schemes into climate models has, for instance, improved global warming projections (Strommen et al. 2019). Theoretically, there are many forms of stochastic processes with memory widely known as anomalous diffusion as exemplified by fractional Brownian motion Biagini et al. 2008;Mishura 2008;Sithi and Lim 1995). The stochastic process with memory used in this paper belongs to a larger class of non-Markovian white noise processes which have been successfully applied recently to investigate other systems such as the ageing of fibrin (Aure et al. 2019), the DNA distribution in genomes (Violanda et al. 2019), diffusion coefficient values for proteins of varying lengths (Barredo et al. 2018), and cyclone track fluctuations (Bernido et al. 2014), among others. We use an analytical stochastic framework with memory Carpio-Bernido 2012, 2015) which allows direct comparison between analytical and empirical results for the mean square deviation (MSD) of observables. Once the theoretical MSD has been matched with the empirically generated MSD of the fluctuating values of the observables, we obtain an explicit analytical probability density function (PDF) describing the GBR degradation, as well as SST and atmospheric CO 2 levels. The PDF has a memory parameter and characteristic frequency with distinguishing values that define the memory behavior and comparative temporal evolution of the GBR degradation, elevated SST and atmospheric CO 2 levels.
The closed analytical expression for the PDF, which satisfies a modified diffusion equation with time-varying diffusion coefficient allows determination of various physical quantities of interest. An added insight into the temporal dynamics of the examined systems is obtained by looking at the short-time behavior of the time series data. Specifically, the fluctuating values of the dataset correspond theoretically to an anomalous diffusion where consecutive values are positively correlated. Short-time behavior of the analytical PDF and MSD reveals that SST and atmospheric CO 2 levels are characterized by fluctuating data point values in the superdiffusive regime, while the GBR degradation exhibits hyperballistic fluctuations.

Observable fluctuations in empirical datasets
We investigate three sets of empirical time series: (a) the mean percentage coral cover of the GBR from July 1985 to June 2011 obtained from the archives of the Australian Institute of Marine Science (AIMS) accessible through the data portal data.gov.au (AIMS 2016); (b) Sea surface temperatures (https:// www. ncdc. noaa. gov/ telec onnec tions/ enso/ indic ators/ sst. php) from the climate monitoring division of NOAA for monthly data from January 1982 to October 2017; and (c) the Keeling curve for atmospheric CO 2 levels measured by the Mauna Loa Observatory of the National Oceanic Atmospheric Administration (NOAA, www. esrl. noaa. gov) in the island of Hawaii. Graphical representations for these three datasets are shown in Figs. 1, 2, and 3. Figure 1 shows the mean percentage hard coral cover of the GBR from July 1985 to June 2011 obtained from the archives of the AIMS. Starting in the early 1980's, these data were gathered using manta tows to survey large areas of reefs and calculating the percentage of the perimeter of each reef that is covered with living and dead corals.
For sea surface temperatures, there are four sets of data available depending on the region in the equatorial Pacific Ocean: Niño 1 + 2, Niño 3, Niño 3.4, and Niño 4. Specifically, we investigate the data closest to the Great Barrier Reef (18.2871° S, 147.6992° E) which is Niño 4 (5°N-5°S, 170°W-160°E) (https:// www. ncdc. noaa. gov/ telec onnec tions/ Our study of the rising atmospheric CO 2 level, on the other hand, considers the Keeling curve for measurements taken since 1958 (www. esrl. noaa. gov) at the Mauna Loa Observatory. In particular, we take the weekly mean CO 2 data from May 1974 to August 2017, shown as open circles in Fig. 3. Whereas fluctuation is quite straightforward for the GBR, and relatively random for sea surface temperatures (Fig. 2), the atmospheric CO 2 levels exhibit a highly oscillatory data and certain periodicity. There are daily CO 2 level oscillations between night and day, as well as seasonal oscillations between Spring and Fall. For the weekly mean CO 2 levels from May 1974 to August 2017, we put a linear fit (solid line), or quadratic fit (dash-dot) for the empirical data (open circles) as shown in Fig. 3.
One can plot fluctuations in the data of Fig. 3 by subtracting each data point from the linear fit (solid line). Likewise, for comparison, we can do the same using the quadratic fit in Fig. 3. The deviations from the linear and quadratic fits are shown in Fig. 4.
From the empirical data shown in Figs. 1, 2, and 4, the mean square displacement (MSD) of the fluctuating observables are computed as a function of time. If the value of an observable is represented by x(i) where i is a specific time, the MSD is computed as, with N the total number of data points. Here, Δt < N , is the time window, or lag time, which is an interval between two data points. Equation (1) gives an MSD value corresponding to a given Δt and to generate a MSD versus time plot, the Δt is increased up to a desired value. Note that as Δt increases, the number of sample deviations to be averaged in the dataset decreases. To avoid a lack of statistics, an upper limit for the value of Δt is imposed. This empirical MSD versus time plot could then be matched with a theoretical MSD of the stochastic analytical framework.

Stochastic framework with memory
The dataset we are using exhibits fluctuating values of the measured observables. Mathematically, we can capture the behavior of these fluctuations using a modulated Brownian motion B(t). Designating a fluctuating variable as x(t) which starts at x 0 at time, t = 0 , and ends at, t = T , we can explicitly write, x(t) = x 0 + fluctuations , by parametrizing it as, In Eq. (2), the factor (T − t) ( −1)∕2 endows the system with memory as it evolves from time, t = 0 , to a final time, t = T = ∕2 . The memory behavior of the system changes depending on the value of the memory parameter . On the other hand, the factor t ( −1)∕2 √ cos( t) in Eq.
(2) modulates the Brownian motion B(t) where is a characteristic frequency of the fluctuating measured observables. As shown in the Appendix, if we fix the endpoint of x(t) to be at x T at final time t = T , Eq. (2) leads to a probability density function (PDF) of the form (see Table 3.1, page 25, of Bernido and Carpio-Bernido 2015), In Eq. (3), Γ( ) is the Gamma function and J (z) is the Bessel function of the first kind. With this PDF one can also evaluate the mean square displacement (MSD) for x(t) , i.e., x, T;x 0 , 0 � dx , with P x, T;x 0 , 0 given by Eq. (3). This yields the result (Bernido and Carpio-Bernido 2015), The values for the parameters and would depend on the system being studied.
We observe that the probability density function, Eq. (3), can be written in the form, with a memory-dependent MSD, Eq. (4). This PDF gives us the flexibility of determining an initial and final point for the evolution of a system even for large times. Moreover, Eq. (5) obeys the modified diffusion equation (Barredo et al. 2018;Bernido and Carpio-Bernido 2015), Equation (6) exhibits a time-dependent diffusion coefficient of the form, (MSD)∕ t , which interestingly also occurs in other physical phenomena (Lee et al. 2020;Reynaud 2017;Li et al. 2018;Jeon et al. 2014;Aure et al. 2019;Barredo et al. 2018). We now proceed to apply this stochastic process with memory to acquire insights into the empirical datasets for the GBR degradation, sea surface temperatures near the GBR, and the rise of atmospheric CO 2 levels.

Great Barrier Reef: empirical and theoretical mean square displacements
To check whether the theoretical MSD, Eq. (4), indeed describes fluctuations for the percent of hard coral cover The normalization N effectively raises the MSD graph with its shape unaltered, and T C shifts the whole plot to the right so that empirical and theoretical plots are superimposed for easier comparison. 1 3

Sea surface temperatures: empirical and theoretical mean square displacements
We now consider fluctuating values of sea surface temperatures shown in Fig. 2. Evaluation of the mean square deviation yields an MSD shown in Fig. 6. For SST, a match for the empirical MSD is obtained using the same stochastic behavior as the GBR degradation, i.e., Eqs.
(2) to (4). In particular, empirical and theoretical MSD plots match for parameter values, = 1.18 and = 0.19 , as shown in Fig. 7. Again, for ease in comparing empirical and theoretical plots, we lower the analytical MSD graph without altering its shape using a normalization N (MSD), where N = 0.11.

Atmospheric CO 2 levels: empirical and theoretical mean square displacements
Given the empirically derived fluctuations shown in Fig. 4 for atmospheric CO 2 levels, one can plot the corresponding mean square displacements as shown in Fig. 8. The cyclic nature of atmospheric CO 2 levels clearly manifests in its log-log plot of MSD. For CO 2 levels, empirical and theoretical MSD plots match for = 1.25 and = 0.07 . For ease in comparing empirical and theoretical plots, we lower the analytical MSD graph without changing its shape with the normalization N (MSD), where N = 0.4. The comparison is shown in Fig. 9.

Discussions
To help understand the similarities and interrelated dynamics of atmospheric CO 2 levels, SST, and declining hard coral cover of the GBR we used a stochastic process with memory as an analytical framework. A good match between  empirical and theoretical MSD's was attained for each phenomenon indicating that all three systems are described by the same class of non-Markovian stochastic processes described by Eqs.
(2)-(4). A comparative insight into the complex dynamics of the GBR degradation, SST, and atmospheric CO 2 levels can be seen by summarizing values of the memory parameter μ and characteristic frequency as shown in Table 1. For both the memory parameter and characteristic frequency , the GBR degradation registers the highest values. Although the three systems are influenced by the same memory behavior of the form, (T − t) ( −1)∕2 , the effect of the memory function differs with values of = 4.64 for the GBR, = 1.25 for CO 2 levels, and = 1.18 for SST. Given the same values for an interval (T − t) , the impact of this memory range is greater for larger values of µ > 1. Also, as seen in Eq. (2) the stochastic process loses its dependence on the memory function, (T − t) ( −1)∕2 , at the critical value of = 1.
For short-time intervals, all three systems exhibit a memory behavior different in form from Eqs. (2) to (4). Note that for T ≪ 1 , we can use the asymptotic forms, cos( T∕2) ≈ 1 , and J −(1∕2) ( T∕2) ≅ 1∕Γ( + (1∕2)) ( T∕4) −(1∕2) , which reduce Eq. (4) to a mean square deviation exhibiting a power law, i.e., where, = 2 − 1 , and c = √ 4 Γ( )∕4 Γ( + (1∕2)) . The MSD, Eq. (7), is typical of anomalous diffusion with exponent (Safdari et al. 2017;Wang et al. 2020;Metzler et al. 2014). Table 1 gives us the values of the memory parameter which yield, = 8.28 for the GBR, = 1.50 for CO 2 levels, and = 1.36 for SST. The value of has a physical implication on the behavior of the fluctuating values of the observables. In particular, 0 < < 1 , gives subdiffusive behavior whereas, 1 < < 2 , is superdiffusive Biagini et al. 2008). The case = 2 is the ballistic domain, while 2 < is referred to as hyperdiffusive or hyperballistic growth of the MSD that transpires, for example, in particle diffusion in turbulent flows or hypertransport in optical systems (Safdari et al. 2017;Wang et al. 2020;Peccianti and Morandotti 2012;Levi et al. 2012;Siegle et al. 2010;Yao et al. 1995). For, 1 < , there is long range memory behavior where a fluctuating variable quickly departs from its point of origin as time evolves in comparison with ordinary Brownian motion. Given these values of the anomalous exponent for T ≪ 1 , the fluctuating data points for CO 2 levels and SST belong to the superdiffusive regime, while the values for the GBR is hyperballistic in its growth of the MSD. Note that, although the short-time MSD, Eq. (7), is anomalous in nature, the stochastic process used here is Gaussian. This should be differentiated from (7) MSD ≈ cT , the anomalous but non-Gaussian phenomenon which may involve, for example, the generalized grey Brownian motion (Ghosh et al. 2016;Sposini 2018;Bock 2020).
We illustrate the different fluctuations or diffusion regimes in Fig. 10 by plotting the probability density function with an MSD given by Eq. (7)   x − x 0 are narrowly concentrated about the origin indicating that the value of x hardly fluctuates from an initial value x 0 as shown by its probability density function in Fig. 10. All plots in Fig. 10 were done at T = 0.2 . Note that the critical value, = 1 , or equivalently = 1 in Eq.
(2), corresponds to no memory or ordinary Brownian motion. Figure 11 shows the MSD plots of Eq. (7) for anomalous diffusion ( = 1.5 or = 1.25 ) given by the straight red line, and ordinary Brownian motion for the case with no memory ( = 1 or = 1 ) as shown by the straight dashed line. For comparison, we also graph in Fig. 11 the MSD, Eq. (4), for CO 2 levels with = 1.25 and = 0.07 which coincides with the anomalous diffusion for smaller values of T. Clearly, the memoryless Brownian motion (dashed line) and the power-law form of the MSD, Eq. (7), for anomalous diffusion (solid red line) are not able to capture the dip in MSD for longer times T exhibited by both the empirical MSD and theoretical MSD, Eq. (4), as shown for example in Fig. 9. The MSD plots for the GBR, sea surface temperatures, and CO 2 levels given by Figs. 5, 7, and 9, respectively, all experience a dip at longer time T which may imply a subdiffusive behavior. However, other intervening factors could allow its transition back to a superdiffusive regime. For possible oscillations between diffusion regimes, the GBR could be more prone to this given its high value for the characteristic frequency .
The evolution in time for fluctuations corresponding to each phenomenon can be compared by plotting Eq. (3) for PDF versus time as illustrated in Fig. 12. Here we set, Δx = x T − x 0 = 2 , in Eq. (3) and use the values for µ and ν in Table 1. Figure 12 shows that the PDF for all three phenomena reaches its peak roughly around the same time. The SST and atmospheric CO 2 levels essentially have the same probability distribution that appears to slowly diverge at longer times. The PDF of the GBR, however, drastically drops right after reaching its peak. This scenario implies that for a difference of 2 percent between initial and final hard coral cover, i.e., Δx = 2 , the probability that this difference persists is close to zero after a certain time. On the other hand a difference in fluctuation of 2 °C for SST and 2 ppm for CO 2 levels are possible and sustained even for longer times as depicted in Fig. 12.

Conclusion
The unified stochastic framework with memory gives us a snapshot through a common lens of three seemingly disparate but interrelated complex systems: the GBR, SST, and atmospheric CO 2 . The inclusive framework allows for a comparison of the patterns exhibited by the datasets generated by each system. The translation of the behavior of datasets into an analytical stochastic setting shows that all three systems possess long term memory properties. Equation (7) tells us that for short-time intervals, T ≪ 1 , the datasets for all three systems reveal a mean square deviation of the measured observables growing in time as, MSD ≈ cT , where = 2 − 1 , and c depends on the memory parameter (see Table 1). Each dataset examined exhibits anomalous diffusion where successive fluctuating values of measured observables are positively correlated. Specifically, for atmospheric CO 2 and SST, we have superdiffusive values of the exponent where, 1 < < 2 , while for the GBR degradation the value is, 2 < , which falls in the hyperballistic domain. Stochastic acceleration or hyperballistic diffusion has been shown to result from a dynamically evolving disorder as demonstrated experimentally and numerically using an optical setting (Levi et al. 2012). The lesser values of the anomalous exponent for atmospheric CO 2 ( = 1.50 ) and SST ( = 1.36 ), as compared to a large value of = 8.28 for the GBR, may be attributed to the fact that atmospheric CO 2 and the ocean are much larger systems than the GBR and, therefore, more capable of absorbing ecological disturbances. Their relative values, however, quantify the observation that small fluctuations in atmospheric CO 2 and ocean temperatures could have a dramatic effect on the GBR. At finite times T , the MSD for the three systems transition to the full expression given by Eq. (4).
From the MSD, Eq. (4), characterized uniquely by the memory parameter and characteristic frequency for each system (Table 1), we also obtained an explicit form of the probability density function. The explicitly derived probability density functions and the modified diffusion equation, Eq. (6), for degradation of the Great Barrier Reef, sea surface temperatures, and atmospheric CO 2 levels provide a novel perspective and new insight on the temporal interrelated dynamics of these systems and a mathematical handle on effects of collective ecological memory.

Appendix: Deriving the probability density function
Consider a fluctuating variable x(t) which starts at x 0 at time t = 0 . Using the Brownian motion B(t) , we parametrize fluctuations up to some final time t = T given by Eq.
(2) as, where, (T − t) ( −1)∕2 , is a memory function, with t ( −1)∕2 √ cos( t) modulating ordinary Brownian motion B(t) . To evaluate the probability density function for this stochastic process, we express dB(t) as dB(t) = (t)dt , where (t) is the derivative of the Brownian motion (Hida et al. 1993), i.e., (t) = dB(t)∕dt . Note that, Eq. (8) fixes the initial point at x 0 but the fluctuating variable can end anywhere at some final time t = T . We can fix the endpoint of x(t) at time T to be, x(T) = x T , using the Donsker delta function, x(T) − x T . This delta function is always zero for all paths described by Eq. (8), except those which end at the chosen endpoint x T . The expectation value for all contributing paths ending at x T is an integral over the white noise variable (t), where d ( ) is the Gaussian white noise measure (Hida et al. 1993;Hida and Streit 2017;Kuo 1996;Obata 1994), and P x T , T;x 0 , 0 is the probability density function. To evaluate Eq. (9), we write the Donsker delta function in terms of its Fourier representation, (x) = (1∕2 )∫ +∞ −∞ e ikx dk , to get, Using Eq. (8) for x(T) we have, Integration over d ( ) can be done using the characteristic functional (Hida et al. 1993), where we let, (t) = k(T − t) ( −1)∕2 t ( −1)∕2 √ cos( t) in Eq. (11). The key to making the model computable in the sense of explicitly evaluating expressions involving the Gaussian white noise variable (t) is the characteristic functional, Eq. (12), which defines the Gaussian white noise measure d ( ) (Hida et al. 1993;Hida and Streit 2017). Moreover, Eq. (12) allows the generation of a wide class of stochastic processes with memory which has been applied to other biological and physical systems (Aure et al. 2019;Violanda et al. 2019;Barredo et al. 2018;Bernido and Carpio-Bernido 2015). With Eq. (12), the probability density function, Eq. (11), becomes, The remaining integral over dk is a Gaussian integral which yields the probability density function (Bernido and Carpio-Bernido 2015), We then obtain Eq. (3) using Eq. (3.768.9) of Gradshteyn and Ryzhik (1994) to integrate over dt in Eq. (14).
Acknowledgements The authors are grateful to Yutaka Shikano for his valuable comment. A. R. E. wishes to thank the Commission on Higher Education and C. B. C. the Department of Science and Technology for funding support. C. C. B. and M. V. C. B. acknowledge the support of PLDT-Smart Foundation.

Conflict of interest
The authors have no conflicts of interest to declare that are relevant to the content of this article. (11) .