Generalization of Fluctuation-Dissipation Theorem to Systems with Absorbing States

Systems that evolve towards a state from which they cannot depart are common in nature. But the fluctuation-dissipation theorem, a fundamental result in statistical mechanics, is mainly restricted to systems near-stationarity. In processes with absorbing states, the total probability decays with time, eventually reaching zero and rendering the predictions from the standard response theory invalid. In this article, we investigate how such processes respond to external perturbations and develop a new theory that extends the framework of the fluctuation-dissipation theorem. We apply our theory to two paradigmatic examples that span vastly different fields - a birth-death process in forest ecosystems and a targeted search on DNA by proteins. These systems can be affected by perturbations which increase their rate of extinction/absorption, even though the average or the variance of population sizes are left unmodified. These effects, which are not captured by the standard response theory, are exactly predicted by our framework. Our theoretical approach is general and applicable to any system with absorbing states. It can unveil important features of the path to extinction masked by standard approaches.


I. INTRODUCTION
Statistical physics has been an invaluable bridge between microscopic dynamics and macroscopic observables.The Einstein-Smoluchowski equation, the regression hypothesis by Onsager [1][2][3], and linear response theory by Kubo [4,5] are different kinds of fluctuation-dissipation theorems (FDTs).These relations, which were also extended to stochastic systems [6,7], connect the response of an observable to stationary correlations.This allows for the prediction of the effects of perturbation from equilibrium dynamics.Though equilibrium dynamics is a useful description, studies into non-equilibrium systems have gained much traction due to the vast number of examples showing non-equilibrium behaviour.These examples include processes from biology, epidemiology, turbulence, and chemical systems, amongst others [8][9][10][11].
Given the predictive power of these theorems, their violations also become crucial to study.A substantial body of literature exists on the violation of the fluctuation-dissipation theorem in glassy systems [12,13].In non-glassy systems, there have been attempts at extending the validity of FDT to non-equilibrium dynamics using the concepts of asymmetry [14,15], frenesy [16], and local currents [17][18][19].In particular, successes have been obtained at predicting response to perturbation in an arbitrary non-stationary state [20].However, the prevalent and inherent assumption is that there is some form of eventual non-trivial steady state.Instead, the focus of this work is on an entirely different class of non-equilibrium systems with absorbing states/boundaries.
Processes with absorbing states are ubiquitous in nature.Some examples include chemical reactions, epidemics, and population dynamics [21][22][23][24].In these cases, the net outward flux of the system's constituents is positive, and hence the total probability distribution decays with time, rendering the steady-state trivial.However, they are also frequently found in a quasi-stationary regime for a long time before reaching extinction.There exists extensive mathematical literature analyzing the properties of quasi-stationary systems [24][25][26].The details of time taken to reach extinction, the existence of quasi-stationary distribution [27,28], simulation methods [29], etc have been studied with applications to cellular automata [30], birth-death processes [31], Brownian motion [32], contact process [33], and many others.
Despite the vast literature, there have been no forays through the lens of statistical physics into looking at the change in FDT in the presence of an absorbing state.There is a gap in our understanding of the effects of perturbation in decaying processes.Linear response theory in "quasi-stationary" states has been previously considered [34,35], but the states referred to are metastable states in which systems spend a long time before reaching equilibrium, rather than processes reaching extinction.In general, fields of quasi-stationary systems and non-equilibrium FDT are both crucial and significant, but they have not been connected.
In the following work, we build the bridge between these two fields by considering the modifications needed by linear response theory to accommodate the presence of absorbing states.We consider the method of conditioning observable averages over trajectories that have not yet hit the absorbing state, calling this operation conditioning to survival.With this framework, we find that a generalized FDT arises by introducing a new observable that accounts for the survival of trajectories.This theorem is a true generalization and cannot be obtained by any ad hoc replacement of the stationary distribution in the standard case.
To demonstrate the scope of the approach, we illustrate applications in two different and wide-ranging contexts (ecological and biochemical), both having important practical implications.
One of the standard techniques used for stochastic modelling of forests is the birth-and-death process [36,37].This one-step Markov process, when considered with constant per capita birth rate (b), death rate (d) and immigration, results in the well-known Fisher log-series stationary distribution [38] which gives the relative species abundances (RSA) and a biodiversity parameter in an ecosystem.Under the absence of immigration, including the zero population absorbing state (and considering b < d), we show the validity and importance of our theory, comparing it to the equivalent of the standard case, where the response is significantly underestimated.We then proceed to show that, while common observables like population size average and variance, and even the observed RSA can remain constant under perturbation, the extinction probability for a species can irreversibly be changed.
On the other scale, we pick the example of targeted search by proteins on DNA.This process is modelled using hybrid diffusion models where the protein moves along the DNA searching for a specific set of base pairs (1D diffusion) but can also detach temporarily and move in free space (3D diffusion) [39,40] leading to speeding up of the target finding rate [41][42][43].In this example, we consider a discrete state version of 1D diffusion of the protein, with a possible detachment to 3D space and reattachment to the DNA.This kind of discrete model has been used previously to explain different properties of the process, especially optimal target-finding speed [44,45].Here, we predict the distance of the protein from the target under changing temperature conditions and proceed to show how the distribution of absorption time skews towards faster absorption with a sudden perturbation in the detachment rate.
In both examples, we compare analytical predictions to numerical simulations of the examples and find excellent agreement between the two.

II. GENERALIZATION OF FLUCTUATION-DISSIPATION THEOREM
We begin with a recap of the standard FDT that we will later use to compare the generalized version.Consider a stochastic system comprised of a finite set of discrete states, represented by Ω.Let P x (t) be the probability of finding the system in state x ∈ Ω at time t.The evolution of the probability P x (t) is governed by the Master equation where W xy is the rate of transition from state y to state x.By defining a matrix with elements H xy ≡ W xy − δ xy z W zy , the Master equation can be written as and hence has the solution P (t) = P (0) e H t [24].Without any absorbing boundaries, at large times, if the process is irreducible and positive recurrent [46], a stationary state is reached, where the probability of finding the system in state x is P st x .A perturbation is described by H + δH(t), with δH(t) = F (t)δH where F (t) is the time component of δH.If the transition rate matrix H depends on a certain set of parameters f ≡ (f 1 , f 2 , . . .), the perturbation changes at least one of them, resulting in f → f + ∆f (t).Then, δH(t) = (∂ f H) ∆f (t) at leading order and the response of any observable A is given by the derivation from its average stationary value according to [6] x,y,z where ⟨• • • ⟩ represents the average over the probability distribution, R A,f (t) is the response function which determines the response of A to perturbation f + ∆f (t), Θ(t) is the Heaviside step function, and ϕ x ≡ ln P st x is the potential characterising the stationary solution (See Supplementary Material Section 1 for detailed derivation).This equation connects the fluctuation of two variables at equilibrium to the dissipation back to equilibrium after a perturbation by expressing the response function as a time derivative of a two-time correlation function in unperturbed dynamics.Hence, such a relation is termed the fluctuation-dissipation theorem.
In case the process has absorbing states, they are characterised by outgoing transition rates being zero in the master equation, i.e., W x,y = 0 ∀x ∈ Ω , y ∈ ∂Ω, where ∂Ω the set of absorbing states.We are interested in how the observables can be measured in the interior Ω • = Ω − ∂Ω.We collapse all the absorbing states into a single representative state o to avoid degeneracy.We are then able to write the master equation in the same form as Eq. ( 2), but with redefined transition rates on Ω ′ = Ω • ∪ {o}.The right eigenvectors and eigenvalues of H will be denoted as φ n and λ n .The entries of the first column of H are all zeros; hence the first eigenvalue is λ 0 = 0 and the corresponding first right eigenvector is φ 0,x = δ x,o .This corresponds to the long-time limit of the probability distribution, which is zero everywhere on the interior.We consider the matrix W to be an irreducible matrix, which means there always exists a non-trivial path from any interior state to any other interior state.From the Perron Frobenius Theorem [47], the eigenvalue spectrum is where ℜ(•) denotes the real part of the argument.The second eigenvalue λ 1 is real and the corresponding left and right eigenvectors are real with positive components in the interior (See Supplementary Material Section 2.1.1 for proof and details about the heirarchy.See also [48,49]).
Then, the propagator giving the probability of being in the state x at time t starting from y at time t ′ can be written in terms of the eigenvectors and eigenvalues of H, i.e., P (x, t|x 0 , 0) = n c n,x0 φ n,x e −λnt where c n,x0 are the constants depending on the initial condition.For times t such that t(ℜ(λ 2 ) − λ 1 ) >> 1, the probability distribution is approximated by P lt x (t) = φ 0x + c 1 φ 1x e −λ1t .To ensure that the system has not reached the boundary yet, we introduce an observable that characterizes the survival, namely, χ x = 1 ∀x ∈ Ω • and χ o = 0.The conditional average of an observable A on the interior (i.e., A o = 0) reaches a stationary value as t → ∞, i.e., where we have used the eigenvector expansion for the propagator and φ 1x is the eigenvector associated with λ 1 (See Supplementary Section 2.1.2for more details).It is clear from this definition that we are conditioning the observable to survive on the interior.The superscript lt is used to denote the conditional observable average at large times.Similarly, we write the conditional correlation between two observables A and B at times t and t ′ , with t > t ′ , which at sufficiently large starting times gives The averages in Eq. ( 7) and Eq. ( 8) are equal to averaging the observables directly over P lt x (t) which is self-consistent since both are considered at long times.
A perturbation in the parameters can be represented by a perturbation in the generator H → H + δH(t), leading to a change in the long-time probability distribution, P (t) = P lt (t) + δP (t).The solution from the first order differential equation for δP is, The change in an observable conditioned to survival is − ⟨A||χ⟩ lt δ⟨χ⟩ t ⟨χ⟩ lt t + higher order terms (10) which can be obtained from fraction differentiation.We have used the superscript lt to denote that we are starting from the observable average approximately given by Eq. (7).Assuming that the perturbation can be factorised, i.e., δH(s) = V δF (s) and B x ≡ 1 φ1x y V x,y φ 1y , the response function is a correlation function.Considering a perturbation affecting a system parameter, we obtain the generalized fluctuation-dissipation theorem similar in form to a standard FDT (see SI Section 2.2 for detailed derivation): where φ is now the potential associated with the second right eigenvector, rather than the stationary state, i.e., φ 1,x = e φx .The difference between the standard case is immediately apparent from the appearance of the second term which does not exist in Eq. ( 5).In the standard case, this term can be shown to be zero.Consequently, while the standard case and the absorbing boundary case cannot be directly compared, the first term corresponds to the standard FDT with the quasi-stationary distribution used in the place of steady-state distribution.The additional term arises from the fact that the standard FDT considers all states but we need to consider only the survived trajectories, and hence, only the interior states for cases with absorbing states.This term is important as it ensures that the average is performed over trajectories that survive at least until time t which is the time of our observation.
The response function given in Eq. ( 12) is valid for all state-dependent variables, such as mean, variance, skewness, etc.But since Eq. ( 9) is valid for all absorbing processes, we can use it to compute the change in probability distribution dependent quantities such as survival probability and first passage time distribution.The first passage time is defined as the time taken for a system to reach the absorbing state and can range from zero to infinity.Even in the standard theory, an absorbing boundary equivalent state is needed to define these quantities.These state-independent quantities are also significant in many practical examples.
The survival probability S(t) and the first passage time distribution P(t) are given by Usually, the mean first passage is computed for ease of calculation, but having the full solution for the perturbation and the response function, we can compute the change in the entire first passage time distribution.The resulting response functions for the survival probability and the first passage time distribution due to a perturbation in the parameter f are Equations ( 12), (15), and ( 16) are our main results.Although they are general and strictly valid for all finite and discrete systems, they can be extended to infinite systems that obey the same eigenvalue spectrum as in Eq. ( 6).We proceed to show certain applications of these response functions in systems of interest.

III. APPLICATIONS A. Birth-Death Process
The birth-death process is a common starting point for modelling stochastic population dynamics in ecological communities.A random fluctuation may lead a species to reach zero population, from where it will never recover (without immigration).To investigate how populations approach extinction, we consider a model where the effective (per capita) rates of reproduction and death of an individual are constant.Thus, the rates of the master equation are given by W n+1,n = b n = bn (birth event), and W n−1,n = d n = dn (death event).Here the population size of a species is given by the discrete variable n = 0, 1, 2, . . . .Under the neutrality assumption [38,50], the probability distribution of n gives the distribution of the population of all species in the system.When b < d the state n = 0 is an absorbing boundary.The eventual probability that all species' populations reach zero is one since there are no possible birth events when a species has gone extinct.
A key point is that the birth-death process is an infinite state system, but still satisfies Eq (6), because the eigenvalues are discrete and are given by λ k = (d − b)k for k ≥ 1.Hence, we proceed to use the derived response theory and FDT since the necessary hierarchy is still satisfied.With b −1 = 0, the Master equation is given by Ṗi The second eigenvector and the eigenvalue can be obtained by the generating function [22] G(z, t) At large times A(z, t) is small enough to use the approximation 1/(1 − x) y ≈ 1 + yx.Since the coefficient of the slowest exponential in the long-time state gives the second eigenvector, we obtain φ 1n = ν(1 − ν) n−1 .The time-dependent solution to the Master equation is calculated by using the Meixner Polynomials [36].(See Supplementary Material Section 3 for the full form of the time dependent solution).
Figure 1a shows the comparison between results from simulation, the new theory given by Eq. ( 12), and the standard-equivalent case given by only the first term of Eq. ( 12).The first term of Eq. ( 12) is similar to Eq. ( 5), but with the quasi-stationary distribution in place of the stationary distribution, corresponding to a direct substitution in standard theory.We consider the average population squared, i.e., ⟨n 2 ⟩, since this shows a significant deviation from standard theory (for details about simulations, see Appendix).It is immediately apparent that the standardequivalent case significantly underestimates the response of the system while the predictions from the new theory and results from simulations match very well.This indicates that the direct replacement of the distribution is a wrong approach and could result in an incorrect estimate of the response.FIG. 1. Simulated and predicted responses in the birth-and-death process with absorbing boundaries: a) Comparison between modified response theory Eq. ( 12), standard-equivalent case (similar to Eq. ( 5)), and simulations of the system.The observable considered is the average population size squared (An = n 2 ).The plots show the deviation of survived observable from the unperturbed value.The green line indicates the contribution from the first term in Eq. ( 12).The black line is the prediction using both terms, and the majenta line is the deviation observed in simulations.The perturbation strength is ∆b = b/10.b) Increase in extinction probability of a species due to sustained perturbation for a period of t = 1 year in both birth and death rates.Inset shows the average survived population (from simulations) against time in both perturbed and unperturbed cases.An important quantity for ecosystems is the survival probability of species.An increase or decrease in this quantity could mean the difference between extinction and persistence.To study this, we consider a sustained period of perturbation in both rates, reasoning that any increase in death rate caused by natural or man-made causes results in more available resources, thereby also increasing the birth rate.After computing the response function for survival probability, it is possible to compute the extinction probability which is given by 1 − ⟨χ⟩ t .
Figure 1b shows the extinction probability for a perturbation in both the birth rate and the death rate that lasts for a period of t = 1 year while comparing it to the extinction probability in the unperturbed case.The inset shows the change in average survived population size during and after the perturbation.It is seen that the average population size is unaffected by this perturbation.While not plotted, the quasi-stationary distributions in both cases remain the same, despite the perturbation, which indicates that the observed relative species abundance distribution can appear to be constant in time.The critical point of note is that the extinction probability of the species immediately increases, especially during the period of perturbation, and remains larger than the unperturbed case even after the perturbation has ended.The increase in extinction rate caused by the changing population distribution in both cases is exactly compensated by the changing survival probability, which leads to the constant quasi-stationary distributions, but the increase in death rate is not aptly compensated by the birth rate increase and hence, more species could go extinct.This means that even though the observed average population size and the abundance distribution of the system stay constant, the total number of species in the system decreases faster than usual, and this effect persists beyond the perturbed duration.
While the constant population size and the relative species abundance is ultimately an effect of the particular choice of same perturbation strength for the birth and death rate, the observed phenomenon of increased extinction probability demonstrates the importance of investigating ecological systems using more comprehensive statistical tools.

B. Targeted Search on DNA by Proteins
We move from the ecological time and length scales to the biochemical scale involving proteins and DNA.Proteins search for a specific set of base pairs on a DNA to which they have a strong binding affinity, thereby forming an absorbing site in the dynamics.We consider a discrete state model of the targeted search for the binding site, which has already been described in the literature [44,45].The model consists of L − 1 non-specific sites and one specific binding site on the DNA (called 'target').The protein slides along the DNA with a constant diffusion rate u in either direction, until it reaches the target site.At every non-specific site, the protein can get detached from the DNA strand with the rate k of f (which determines the average sliding length), and this results in the protein being in free space.Since diffusion in free space is much faster than diffusion in 1D, we consider free space as a single site.From this site, the protein can attach itself to any spot of the DNA strand with the rate k on .Due to faster diffusion in 3D, all sites are reachable with equal probability, and hence, k on is the same for all sites.
Under this framework, we consider rates estimated from observed experimental data (see Appendix), and extract the eigenvalues and eigenvectors of the transition matrix to perform the computations.For all the simulations of the system, the protein starts from free space and moves according to the transitions described above.Eventually, every simulation will end when the protein reaches the target.But, at any given time, a fraction of them will still have the protein not bound to the target site (which we shall term survived realizations, equivalent to survived species in the birth-death process).The long-time probability of finding the protein at a particular site in these survived realizations is equal to the quasi-stationary eigenvector of the transition matrix.
Binding and detachment of the protein happen through chemical reactions, whose rates can be approximated by the Arrhenius law, given by rate = const × e −Erate/RT , where E rate > 0 is the activation energy of the reaction, R is the gas constant, and T is the temperature of the system.Assuming the Arrhenius law for each of the rates, we consider a periodic temperature change of 5 • C (from 300K to 305K).This results in a perturbation of the rates, which changes the distance of the protein from the target.We compute this change through the root mean squared (rms) distance of the protein from the target.Figure 2a shows that even though the observed pattern of change in distance is non-trivial, the new response theory is able to predict it quite well.The specific pattern of whether the protein moves away from the target or towards it depends on the activation energies (see Appendix).In a chosen experimental setting, the activation energies can be first computed and the theory can then be used to predict the change in distance due to temperature change exactly using Eq. ( 12) and Eq. ( 3).
The time to extinction is an essential random variable in absorbing processes.The average of this distribution, also called the mean first passage time, provides important information about the protein reaching the target [45].Figure 2b shows the change in the absorption time distribution from simulations and from theory when there is a delta-perturbation in the rate of detachment, k of f .We see that the probability of absorption at short times immediately increases by a significant proportion.Since a delta-perturbation is considered, at large times, the change in distribution will be minimal, which is seen in the figure where the perturbed and the unperturbed distributions converge at times much larger than the time of perturbation.The flat nature of the unperturbed distribution is due to it being an exponential decay with a characteristic scale much larger than the time which is plotted.The matching between simulations and the theory opens new possibilities in ensemble experiments to control mean first passage time through perturbations.

IV. CONCLUSIONS AND DISCUSSION
Our results apply to a broad class of systems that were previously intractable from the lens of response theory.Conditioning to survival is a common tool in mathematical literature of quasi-stationary processes.The generalized fluctuation-dissipation theorem provides insights into how this conditioning affects the response to perturbations.Specifically, this new insight is needed because by direct replacement of the stationary distribution with the quasistationary, one significantly underestimates the actual response.Though perturbation of finite linear operators has been studied previously, obtaining the transient dynamics and connections to statistical mechanics and relevant examples are missing to the best of our knowledge.In order to take into account the inherent decay of the system and to predict the response, the generalized Eq. ( 12) is necessary.
The two examples we have chosen also demonstrate the wide-ranging applicability of including absorbing states explicitly into the dynamics.These systems are also important from a practical viewpoint.Loss of diversity due to extinctions is a concerning problem.Our results indicate a need for a more robust understanding of decaying ecosystem models by demonstrating that there could be inherent perturbations in the system that increase the extinction probability while leaving often-used population measures intact.This is especially important at small population numbers where fluctuations could drive some species extinct permanently.
At the other end of the size scale, various models have been constructed to account for the variety of interactions that occur in protein search on DNA.With advances in technology, the ability to experimentally verify these models is an exciting new line of research.A potential use of response theory is to design a perturbation to observe a specific response.Taking into account the effect of absorbing states and designing a perturbation to observe a specific first passage time distribution could lead to many potential biomedical applications.Some examples include the release of specific patterns of concentration of drugs, control of resource capture by chemicals or nanoparticles, mediating the response of microbial systems through controlled manipulation of growth media, etc.
Absorbing processes form a large class of systems which are ubiquitous in nature.Although we have presented the results in the regime of quasi-stationarity, our theory can be easily generalized to systems that have not yet reached quasi-stationarity.This further increases the scope of the work and has the potential to open other doors of investigation, both on the theory and application fronts.

Details of Simulations
For the birth-death process, we start all species with an initial population of n 0 = 100, using the rates mentioned in the figure captions and simulate until the quasi-stationary regime is reached.Then, a perturbation is applied to the birth rate (which lasts for dt = 10 −5 in the case of delta perturbation or sustained for a given period that is mentioned).The response of survived species is then tracked over time.This process is repeated for 10 6 trials to obtain averages.
For the targeted search on DNA, we use the rates mentioned below and simulate the system till time t = 10/λ 1 , where λ 1 is the leading non-zero eigenvalue and is determined by the transition matrix.For the considered observables in the manuscript, the averages are performed over 10 6 realizations of survived trajectories.

Rates and activation energies in targeted DNA Search
We consider the diffusion rate, u = 2 × 10 5 s −1 , the association rate k on = 0.5 × 10 5 s −1 , and the disassociation rate k of f = 1 × 10 5 s −1 , with L = 100 sites and the target at site 60.The rates have been approximately determined from figures of experimental data of Egr-1 protein which binds to a specific 8 base-pair sequence of DNA [45,51] (rates are chosen at about 400mM concentration of Potassium Chloride (KCl)).To make the calculations easier, we rescale time by dividing all the rates by 10 5 to get rid of the factor.
Assuming that the perturbation occurs in one of the parameters of the system, f From ( 8) and ( 10), we can obtain the following Subsequently, the response function in (7) can now be written as Suppose we denote X f = ∂ϕ(t = 0) ∂f , then, ( 12) can be written as where K A,B (t) is the two time correlation function between observables A and B. Equation ( 13) forms the main Fluctuation Dissipation Theorem we wish to pursue.We are considering the perturbation happening to a stationary distribution.Though we have presented the discrete case here, the same calculations hold for continuous state space with the Fokker-Planck operator being analogous to the transition matrix H [1].
2 Discrete System with Absorbing Boundaries

Master Equation and Preliminaries
Let Ω represent the space of all the states.We can denote the boundary of this space by ∂Ω and the interior by Ω • = Ω − ∂Ω.The transition rates between different states shall be represented by W y→x = W xy .The absorbing boundary condition can be represented by setting the rate of transition out of the boundary states to zero.
The master equation for the time evolution of the probability of being in state x at time t is Let a single state 0 represent all the absorbing states.The corresponding entries for this state that would go into the ME would be Note that these also imply

Eigenvales and Eigenvectors
x Ṗx (t) = 0, i.e, the total probability is conserved (when absorbing states are also considered).We can define the operator H with which we have Ṗ (t) = H P (t) Normalizing the left and right eigenvectors of H such that ⟨ψ n |φ m ⟩ = δ m,n , the first left eigenvector is ψ 0 x = 1 ∀x ∈ Ω for the eigenvalue λ 0 = 0.The corresponding right eigenvector is φ 0x = δ x,0 .In general, for the n th right eigenvector and eigenvalue, with ℜ(z) representing the real part of a complex number z, We will always work with an irreducible W , i.e. for any pair of interior nodes, there is always a path of non zero W 's from one to the other.This implies that the boundary state is not visited.In such a system, (19) can be proved using Perron-Frobenius Theorem [2].Note that the inequality between λ 1 and λ 2 is different from the one between the higher eigenvalues.This is due to the spectrum possibly being complex and hence, having complex conjugate eigenvalues which have equal real parts.Due to the Perron-Frobenius theorem, λ 1 will always be real and hence, have no conjugate counterpart.

Eigenvalue Spectrum
Consider an infinitesimal increment in time ϵ.Then, the Master equation can be discretized to the Markov chain (at small ϵ) Consequently, define T ≡ (I + ϵH) where I is the identity matrix.T is a non negative, irreducible matrix of size N × N where N is the number of states.The eigenvalues of T are simply α i = 1 − ϵλ i .Consider the submatrix of T corresponding to the interior, T * .Due to irreducible system, there exists a path from one state to another, and hence, there always exists a k xy > 0 for which the (W xy ) kxy is non zero.Choosing the largest of these (N − 1) 2 positive numbers ensures that T * is primitive, i.e, (T * ) klargest > 0, if ϵ is sufficiently small.
The characteristic function of the eigenvalues of T * is the characteristic function of T sans the (α − 1) factor.Hence, the eigenvalues of T * is equivalent to the eigenvalues of T .Note that T * xy ≥ 0 when x ̸ = y and T * xx = 1 − ϵ z W zx which is bounded allowing us to choose sufficiently small ϵ such that T * xy ≥ 0 ∀x, y.Thus, we can apply the Perron-Frobenius theorem [2] on the matrix T * .Then,the leading real eigenvalue α 1 > 0. For any other eigenvalue α 2 , then, α 1 > |α 2 |.Then, with ℜ(z) and ℑ(z) representing the real and imaginary parts of z, Correlation functions are defined as ⟨A(t)B(t ′ )⟩ = x,y A x P (x, t|y, t ′ ) B y P y (t ′ ).It follows from ( 24) and the Chapman Kolmogorov equation [3] that This has to be true logically too.Since we are calculating the correlation till max{t, t ′ . . .t ′′ }, the only contributing trajectories have to last at least till time max{t, t ′ . . .t ′′ }.Hence, at time ≤ t ′′′ the trajectories have not yet hit the boundary, i.e. χ(t ′′′ ) = 1.But for the contrary case, this is not true since on the l.h.s.we are requiring that trajectories have not hit the boundary at least till time t ′′′ whereas in the r.h.s.we require that they only survive till a time ≥ t ′′ .This is a consequence of having absorbing boundary.In the standard no absorbing boundary case, χ = 1 trivially and x P (x, t|y, t ′ ) = 1.With the absorbing boundary, we are forced to neglect the 0 state and hence the sum over x ̸ = 0 is no longer unity.If the IC is P lt , In this case, a conditional correlation function, which represents the correlation with only surviving trajectories is written as where I is the identity matrix.In this, we define A xy = A x δ xy and B xy = B y δ xy .(Φ 1 ) xy = φ 1x δ xy .Similar to the conditional average, we can start from arbitrary IC and at large times, it becomes equivalent to averaging over the long time state.Hence, in further calculations, we can just use P lt (t) to determine correlations at large times. 1  At large time limit, t ′ → ∞ and denoting ∆t = t − t ′ , the propagator can be approximated as Using this in (28), By defining ⟨B||χ⟩ lt cond ≡ y B y φ 1y ψ 1 y = ⟨B ψ 1 ||χ⟩ lt .Therefore, at large time limit, Since ψ 1 x > 0, we also have ⟨B||χ⟩ lt cond > 0 if B x ≥ 0. If B = χ, ⟨χ||χ⟩ lt cond = y̸ =0 φ 1y ψ 1 y = y φ 1y ψ 1 y = 1 (since ψ 1 0 = 0 thanks to the normalization of left and right eigenvectors).It is important to note that this is a significant difference from the standard treatment where ⟨A(t)B(t ′ )⟩ st = ⟨A⟩ st ⟨B⟩ st at large time 1 In practice in order to calculate conditional averages at large times, as defined above (see for example (28)), we can use the following formula x,y,z Ax P (x, t|y, t ′ ) By P (y, t ′ |z, t 0 ) P initial (z) x,z χx P (x, t|z, t 0 ) P initial (z) (29) and similarly for multiple time correlation functions.The initial condition P initial at time t 0 is arbitrary as far as P initial (0) < 1.But for our purposes, it is easier to choose P initial (0 Similarly we also write δ⟨χ⟩ t ⟨χ⟩ lt t = ∞ 0 ds R χ,V (t, s) δF (s).Ultimately, with the defined B, we have Notice that, thanks to (31) and ⟨χ||χ⟩ lt = 1, we have lim t→∞ RA,V (t) = 0. Now that we have the form of the response function with conditional correlations, we can try and understand the fluctuation dissipation theorems associated with it.

Fluctuation Dissipation Theorem
Similar to the standard case, let us assume that H depends on some parameter f (can also be a set of parameters.For simplicity, we consider one, but the calculations remain the same for multiple parameters).
We can now Taylor expand H around no perturbation in f We drop the subscript ∆f = 0 for visual simplicity.We see that V = ∂ f H(f ).Since (H + λ 1 I)φ 1 = 0, differentiating, Using the above and ( 28) and ( 38 In the standard treatment, the second term involving χ(t)⟨A||χ⟩ in (43), is absent.Its presence is important since, from, the remark at the end of the previous section, we expect that, in general 3 Birth Death Process The eigenvalues for the Birth Death Process are λ i = (d − b)i [4] which can also be deduced from G(z, t).
We also have the solution for the Master Equation [4].A point of note is that we differ in the absorbing state used in [4].While -1 is an absorbing state in [4], we choose 0 as the absorbing state.Therefore, the solution to the Master equation will be same as in [4] but with the index values shifted by 1 unit.We represent P n,n0 (t) as probability of being at n th state at time t starting from state n 0 at time t = 0. We define The summation over k can be written in a compact form as (de bt −be dt )(be bt −de dt ) bd(e bt −e dt ) 2 Γ(n 0 + 1)Γ(n) where 2 F 1 is the hypergeometric function.In terms of the above definition we have: Response Function Since we are dealing with an infinite dimensional operator, we do not know how e Ht acts on ∂ f φ 1 .Hence, we rewrite (43) slightly.A x e Ht) x,y e λ1(t−s) ∂ϕ 1y (s) ∂f φ 1y Since P (t) = e Ht P (0), with P n (0) = δ n,n0 , we have (e Ht ) x,y = P x,y (t).

∂
where we can use master equation for Ṗx,y (t).While closed form expression for R A,f is difficult to calculate, we can compute it approximately using the explicit formula for P x,y (t) from Eq. ( 47) and (46) .
FIG.1.Simulated and predicted responses in the birth-and-death process with absorbing boundaries: a) Comparison between modified response theory Eq. (12), standard-equivalent case (similar to Eq. (5)), and simulations of the system.The observable considered is the average population size squared (An = n 2 ).The plots show the deviation of survived observable from the unperturbed value.The green line indicates the contribution from the first term in Eq. (12).The black line is the prediction using both terms, and the majenta line is the deviation observed in simulations.The perturbation strength is ∆b = b/10.b) Increase in extinction probability of a species due to sustained perturbation for a period of t = 1 year in both birth and death rates.Inset shows the average survived population (from simulations) against time in both perturbed and unperturbed cases.The dotted black line in both plots indicates the time till which perturbation occurs.Solid black line shows the extinction probability in an unperturbed scenario.Dashed lines are predictions, and colored solid lines are observations from simulations with perturbation strengths ∆b = b/10, and ∆d = d/10.Common parameters are b = 0.4 and d = 0.5.

FIG. 2 .
FIG. 2. Prediction of responses in targeted search by proteins on DNA: a) Top panel: Root mean squared distance of the protein from the target DNA site for periodic changes in temperature.The RMS distance is given in terms of base pairs, where 1bp = 340pm The blue dashed line is the prediction and the solid line is observations from simulations of the protein distance.Bottom panel: Change in the rates of the system assuming Arrhenius rate law.Blue regions indicate the normal temperature (300K), and red regions indicate the higher temperature (305K).The colored dashed lines indicate the rates as given in the legend.b) Prediction of deviation of absorption time distribution with a delta-perturbation of k of f in the system.The histogram of absorption time distribution computed through multiple simulations is represented by the bars with the corresponding probability of absorption on the y-axis.Grey bars indicate the unperturbed absorption time distribution, and the orange bars are the deviation upon perturbation.The solid black line is the prediction from the theory.Perturbation strength is ∆k of f = k of f /10.The unperturbed first passage time distribution is exponential with a large characteristic time, hence appearing flat at plotted scales.The effect of the delta-perturbation is observed only at times much shorter than the characteristic timescale.See Appendix for the parameter values used.