Applications of Dynamic-Equilibrium Continuous Markov Stochastic Processes to Elements of Survival Analysis

In this article, we summarize some results on invariant non-homogeneous and dynamic-equilibrium (DE) continuous Markov stochastic processes. Moreover, we discuss a few examples and consider a new application of DE processes to elements of survival analysis. These elements concern the stochastic quadratic-hazard-rate model, for which our work 1) generalizes the reading of its Itô stochastic ordinary differential equation (ISODE) for the hazard-rate-driving independent (HRDI) variables, 2) specifies key properties of the hazard-rate function, and in particular, reveals that the baseline value of the HRDI variables is the expectation of the DE solution of the ISODE, 3) suggests practical settings for obtaining multi-dimensional probability densities necessary for consistent and systematic reconstruction of missing data by Gibbs sampling and 4) further develops the corresponding line of modeling. The resulting advantages are emphasized in connection with the framework of clinical trials of chronic obstructive pulmonary disease (COPD) where we propose the use of an endpoint reflecting the narrowing of airways. This endpoint is based on a fairly compact geometric model that quantifies the course of the obstruction, shows how it is associated with the hazard rate, and clarifies why it is life-threatening. The work also suggests a few directions for future research.


Introduction
Various non-stationary phenomena, which are studied in the natural/life sciences or engineering and described with solutions of deterministic (ordinary or partial) differential equations, develop on the entire time axis.Examples are a steady-state mode of a non-living system or the dynamic-equilibrium (DE) mode of a living system [1] provided that the latter is mature, i.e. considered in a time interval sufficiently far away from the birth of the system.However, the transition from deterministic differential equations to their stochastic counterparts, which are often associated with multi-dimensional non-homogeneous continuous Markov stochastic processes, is not free of problems.In particular, in the most common formulation, any of these processes can be solely defined at the time points to the right from the initial time point.The present work focuses on the versions of the processes that are defined on the entire time axis.These processes are in focus in Section 2. Section 3 considers the stochastic quadratic-hazard-rate (SQHR) model and deals with applications of DE continuous Markov stochastic processes to survival analyses based on this model.Specifications of the SQHR model for problems in chronic obstructive lung disease (COPD) are discussed in Section 4. Section 5 concludes the work and suggests a few directions for future research.

Invariant Non-Homogeneous and Dynamic-Equilibrium Continuous Markov Stochastic Processes
We denote a stochastic process with ( ) ,t χ ξ where ξ ∈ Ξ is a simple event, ( )  denotes time, and Ξ is the space of simple events.We further let M ρ be the set of non-homogeneous continuous Markov stochastic processes on the Euclidean space n  ( 1 n ≥ ), that satisfy the following two properties.
• All processes in M ρ have the same transition probability distribution.
• This distribution is defined at all time points , s t ∈  such that s t < , and has the density, ( ) , , , s x t y ρ . (The fact that the transition density ρ is the same for all processes in the set is emphasized with the subscript in the notation " M ρ ").
As is well known, the transition density ( )  where ( ) δ ⋅ is the n-dimensional Dirac delta-function.is obvious that a process χ is defined at all t ∈  if and only if both marginal-density functions coincide, as is shown in relation (2).This proves the proposition. The above consideration generalizes the notion of an invariant process for the case where the process does not need be stationary.This generalization goes back to [4].Thus, invariant processes are, in general, non-homogeneous.
The example below is probably the simplest example of invariant processes.
Example 1.Consider the case where 1 n = and the stochastic process χ represents time, t, i.e. ( ) ,t t χ ξ = for all t.Obviously, this process is defined and is continuous at all time points and its marginal density at time t is ( ) and is an invariant process.Its invariant probability density is ( ) ( ) Also note that the transition density in (3) satisfies property (1). One of the most important classes of continuous Markov stochastic processes is diffusion stochastic processes (DSPs).A survey of invariant processes of this type can be found in [2] (Chapter 3).However, no criterion for determination of invariant probability densities for non-homogeneous DSPs was known until recently [5].Here is a brief summary.
We need the notations below.• D ρ is the set of all processes in M ρ , each of which 1) is a DSP with drift n-vector ( ) and T ∇ are the gradient and divergence differential operations with respect to the entries of vector x.

• ( )
, f t y is the so-called Fichera drift.Its entries are We also note that, as is well known (e.g. the paragraph above [2] The above notations allow to summarize the aforementioned criterion of [5] as follows. An invariant probability density of the processes in D ρ , is the solution  , 2 H t y ≡ .In this case, the system described by ( 6)-(7) reduces to ( ) Equation ( 9) has under condition (8) the solution where the function ( ) e t is to be determined.In order to do that, one substitutes (11) into (10) The members of manifold (13) These are indicated in [4] without derivation (see also [5], (15).In contrast, expression ( 13) is derived and, thus, generalizes results ( 14) and ( 15) of [4].
Paper [4] also draws attention to the non-uniqueness of the invariant density.
Indeed, the one-dimensional manifold (13) , then the process is completely described with this density.Consequently, there is no need to involve transition density ρ in the analysis of the process.This fact consistently and substantially simplifies the related theoretical and practical studies. It is also worth to notice that, at any t ∈  , the random variable ( ) This density depends on time t as a parameter.It can also depend on other parameters.If these do not depend on t (e.g. are similar parameters o t and o e in ( 13)), then temporal samples of the process χ can be treated by means of common statistical methods (e.g. the maximum-likelihood technique).If the parameters depend on t, then applicable statistical methods need to be indicated or, if necessary, developed.
A discussion on non-homogeneous invariant DSPs can be found in [5].To our knowledge, in the general case of these processes, there is no method for derivation of exact invariant probability densities.One of the techniques that can provide approximate densities is based on the so-called detailed balance (DB) conditions (e.g.[2], (1.12.13) as well as Sections 3.5.2and 3.5.3).If a DSP is stationary and the DB conditions are met, then the DB solution presents the exact invariant density.Otherwise, the DB solution does not exist but the DB approximation can be regarded as a quasi-stationary one for the corresponding solution.The aforementioned method is presented in ([2], Section 3.5.5;see also Section 3.5.4).It provides a simplified DB approximation for the invariant probability density.Theorem 3.2 in ([2], Section 3.5.5)proves the two-sided estimation for this density when both sides are Gaussian densities.
Apparently, the most important example of the invariant processes is the dynamic-equilibrium (DE) ones [1].The DE probability density corresponding to transition density ( ) , , , s x t y ρ is the one determined with limit relation (see [1], ( 10)) where the limit is uniform in ( ) , t y .Relation (16) presumes that the limit function on the right-hand side does not depend on x.This is the very property that enables one to associate this function with an equilibrium.Since the function generally depends on t, the equilibrium is generally dynamic.Relation ( 16) is inspired by the well-known similar limit relation of R.Z.Has'minskii ([6], (9.12) on p. 139) in the case where the equilibrium is time-independent.One can easily see that the DE density determined with ( 16) is an invariant density.In order to realize that, it is sufficient to apply the Chapman-Kolmogorov In connection with the concept of dynamic equilibrium, one can note the following.
A physical equilibrium is independent of t, i.e. static.A DE need not be static.
Thus, it need not be a physical equilibrium.
A non-living system may have physical equilibrium(s).It may also have modes, which can be interpreted as DEs.These modes are known as the steady-state ones.
A living system cannot have static, t-independent equilibrium(s) but, as rule, has a DE.Another example is a mean-field generalization [10] of the classical, Bogolyubov-Born-Green-Kirkwood-Yvon statistical mechanics.This generalization is free from the thermodynamic-limit assumption, and is more compact and flexible than the classical counterpart.

Application of
One of the aforementioned approaches is proposed in [11].It, in particular, presumes that the hazard-rate-driving independent (HRDI) variables (also known as "risk factors"' or "covariates"; e.g.see [12]) are the entries of some vector n y ∈  described with a particular case of ISODE ( 4) where ( ) ( ) ( ) In this equation, the functions a, A, and B are defined on the entire axis  , and t is interpreted as the age of a person whose survival is described with the HRDI vector y.Notice that ISODE (18) is linear in the narrow sense ([3], Section 8.2) and that it was suggested in [11] four years later the above McKean-ISODE of [8].According to [11], the corresponding hazard rate is (see [11], (2)) where ( ) 0 0 t λ ≥ is the baseline hazard rate and ( ) is the so-called "optimal" trajectory of solutions of (18) (which is, however, not endowed with any modeling representation).( 20)

( )
F t being a n n × -matrix is symmetric and positive-definite uniformly in t.
The advantages of model (19) under condition (18) compared to the well-known Cox proportional hazard-rate model are discussed in detail in ( [11], p. 539).
Moreover ( [11], the text below ( 2)), any stochastic solution y of (18) can follow a deterministic trajectory ( ) t ζ .In more precise terms, this behavior means that any solution y converges to ( ) as time tends to infinity. The discussion in ( [11], p. 539) stresses that typical dependences of hazard rates on the HRDI variables are J-or U-shaped.A particular, exponential case of the J-shaped dependences can be exemplified with the Cox hazard model.Along with this, the U-shaped dependences are meaningful for many HRDI variables, for instance, human body temperature, weight, volume, and surface, blood pressure, serum potassium concentration, serum calcidiol (or calcifediol) concentration and other characteristics.Each of them is in a system-relevant bounded interval.
Values in the middle part of the interval correspond to the lowest hazard rates.However, if a value approaches any of the two bounds, the hazard rate rapidly increases.The simplest version of U-shaped hazard rates is quadratic.Also note that both J-and U-shaped dependences are convex.
Still, it may seem that the quadratic expression ( 19) is nothing but a modeling assumption.However, it has a consistent ground.Indeed, it follows from the exact representation, which is valid under rather mild conditions ([13], Section where the term ∇ is described in the second bullet above Equation ( 4) and Notice that the second multiplier in the integrand in (22) is the Hesse matrix.
Comparing (19) and (21), one obtains the relations T , , , These elucidate a number of topics.
it is only applicable to Equation (18) in the case where functions a, A, and B are independent of t, i.e. solutions of this equation are the Ornstein-Uhlenbeck processes (e.g.[3], Section 8.3).However, it is possible to generalize the "mean-reverting" property to the present settings where a, A, and B are t-dependent with the help of well-known results and the above concept of a DE solution.This can be accomplished in the following way.
Let there exist numbers 0 γ > and for all s t < where ( ) and variance matrix

V t C t s B s B s C t s s y e t y e t t y y
Notably, according to (17), the above mentioned convergence in quadratic mean (see the text above ( 27)) implies convergence in distribution.One of the outcomes of the latter convergence is density ( ) (19) is not sufficiently complete.For example, the lack of modeling representations for the "optimal" trajectory ( ) (see the parentheses in (20)), is not resolved in [11] [14] [15].This gap is partly filled with the results discussed in the text on ( 21)-( 26).In addition to that, comparison of the aforementioned convergence in quadratic mean and the convergence noted in Remark 3 shows that ( ) ( ) i.e., the baseline value ( ) t ζ , or the "optimal" trajectory, of solutions of ISODE (18) is provided by the expectation ( ) DE e t of the DE solution (see (27)).Relation (29) agrees with the difference of ( ) from the zero-drift approximation for y as is emphasized in Remark 3. One can also show that the expectation of the second term on the right-hand side of (19) under condition (29) is, in the infinite-time limit, ( ) ( ) ( ) tr ⋅ is the trace of a (square) matrix.

Specifications of the Stochastic Quadratic-Hazard-Rate Model for Problems in Obstructive Lung Diseases
The SQHR model of [11] [14] [15] can be applied to various areas, whereas we Journal of Applied Mathematics and Physics choose to apply it to clinical trials of treatments against obstructive lung diseases (OLDs), i.e. diseases that cause lower airway obstruction (e.g.[17]).These include asthma, bronchiectasis, bronchitis, and COPD (e.g.[18]).Airway obstruction is a blockage of respiration in the airway (e.g.[17]).For example, the citation below is outlined in ( [19], p. 1342 and Figure 1) and explained in ( [20], pp. 448-449 and Figure 7) in more detail.
"Although the measurements of FEV 1 and FEV 1 /FVC (forced expiratory vital capacity, or volume of air expired between full inspiration residual lung volume) provide a reliable way of diagnosing airflow limitation and classifying COPD severity, they cannot separate the precise contribution of either small-airway obstruction or emphysematous destruction to the airflow limitation in individuals with COPD.However, direct measurements of pressures and flows within the lung indicate that the smaller bronchi and bronchioles less than 2 mm in diameter are the major sites of airway obstruction in COPD.Moreover, the reduced expiratory flow that defines COPD results from reduction of the lumen by peribronchiolar fibrosis, thickening of the small-airway walls, and occlusion of the lumen of the small conducting airways by exudate containing mucus [20]." The above histopathological results on the major sites enable one to reveal the physiological and biological meaning of the hazard rate in the case of OLDs (in particular, COPD) in the compact form.This form, following the definition of the hazard rate, includes relations where t is, as in model (18), the age, 1 m ≥ is the (integer) number of the terminal or respiratory (also known as transitional) bronchioles in the lung (e.g.[21], Section II and Figure 3), ( ) is the area of the cross-section of the lumen (non-occluded or occluded) of the kth bronchiole (e.g.[22], Section "Results" and Figure 1), ( ) λ is the kth-bronchiole hazard rate (assumed to be independent of ( ) . Linear ODEs (30) are complemented with expressions for the survival functions of the bronchioles, ( ) k S t , and the survival function of the entire bronchiole system, ( ) In the particular case where all k λ are the same and independent of t, one can readily check that survival function (32) corresponds to the exponentiated exponential distribution with the parameters independent of t.Also note that expression (32) is well known in reliability theory where it is associated with the so-called parallel systems, and the survival functions are termed the reliability functions.One can also show that this expression corresponds to m stochastically independent random variables.They are associated with the component-individual hazard rates by means of (31).
The role of the bronchiole-specific hazard rates (30) is illustrated in the following remark.
Remark 5.In the particular case where, firstly, the inner surface of the kth bronchiole is circular cylindrical of the non-occluded radius ( ) k R t and, se- condly, this surface is covered by the occlusive layer of highly viscous mucus (e.g. [21], Figure 10 and p. 543) of the thickness  holds and one can readily show that a linear ODE in system (30) applied to the above bronchiole reduces to a linear ODE ( ) ( ) ( ) ( ) ( ) The solution of this initial-value problem is the model corresponds to the core measurement results of [23].The range of the related measurement data is contributed with advances in COPD imaging [24].
These outcomes can help to better focus research on further modeling of hazard rate ( ) t λ (see (30)) for OLDs and interpretation of the results in a less arbitrary way.
The quantity is a biophysical characteristic.Along with this, clinical studies often focus on exacerbations of COPD, which are presumably coupled with the occlusion-caused lumen narrowing (e.g.[20], Section "Acute exacerbation") and are usually formulated in terms of patient symptoms and medical signs (e.g.

test results
).In this case, modeling of ( ) t λ can be based on (19) but the HRDI-variable vector, say, v, can include not only variables on the entire axis  but also variables that are non-negative or represent membership functions (MFs) (e.g.[25]).This means that v is in a bounded set n R ⊆  rather than in the entire Euclidean space n  .However, the latter two types of variables can be represented with entries in  of vector y.Examples of these representations, which presume that a scalar * v is an entry of a vector v and a scalar * y is the corresponding entry of a vector y, are the following: where * 0 s > and * p are the scaling and positioning coefficient.In the latter case, * v can also represent a categorical variable (a particular case of an MF one) such as any of the following three variables (e.g.[26]): • scores 0, 1, 40 on the COPD assessment test; • grades 1, 2, 3, and 4 according to the GOLD; • Since vector-function ( ) , q t v and matrix-function ( ) , Q t v in (37) allow for variables of three different types (see (34)-(36)), Equation (37) can hardly be linear in v.In general, models and methods for ISODEs cannot be applied to (37) directly because ISODEs are usually considered in the entire Euclidean space n  rather than in bounded set n R ⊆  (e.g.[3], Section 6).In order to resolve this matter, one can, in Equation (37), pass from vector v to vector The resulting equation is an ISODE for vector y in the form (4). Remark 6.The generalized SQHR model suitable for methods of DSP analysis comprises nonlinear ISODE (4) for vector y of the HRDI-variable representatives (such as terms * y in (34)-( 36)), expressions ( 21)-(25) under specification (29) and the three properties, which are formulated in the text below (25) and, because of ( 23 • a non-autonomous nonlinear ODE system of the second order ([2], (2.3.9), (2.3.10), and Section 2.3.2) with the initial conditions for entries of the expectation vector ( ) e t ([2], (A.6) and (2.3.12)); the unique advantage of this system is that it includes the influence of diffusion matrix (5) on the expectation; • a non-autonomous linear ODE system of the first order ([2], (2.5.8), (2.5.13), (2.5.14), and Section 2.5.2) with the initial conditions for entries of variance matrix ( ) V t ([2], (1.6.14), (1.6.8)).
The DE versions of the expectation and variance are provided by numerical integration of both systems at time intervals, which is sufficiently far away from the initial time point.One can use a marginal probability density for y, which is the DE density and Gaussian, with the corresponding DE expectation and variance.These outcomes are approximate representations.Reference [2] includes other approaches that are more precise, such as the ones based on the Schauder bases of function Banach spaces and differential-quadrature/stochastic-adaptive-interpolation method.

Concluding Remarks
The above analysis of the SQHR model of [11] [14] [15] generalizes the reading of its Itô ISODE for the HRDI variables.Moreover, it specifies key properties of the hazard-rate function.In particular, it reveals that the baseline value of the HRDI variables is the same as the so-called "optimal" trajectory in the SQHR model and is the expectation of the DE solution of the ISODE.The work also suggests practical settings for obtaining multi-dimensional probability densities necessary for consistent and systematic reconstruction of missing data by the Gibbs sampling, and further develops the corresponding line of modeling.The resulting advantages are emphasized in connection with general survival analysis and the statistical framework for clinical trials of new treatments.The present work also proposes a use of endpoints reflecting the narrowing of airways, which is a major component of obstructive lung diseases (including COPD).This endpoint is based on a fairly compact geometric model that quantifies the course of the obstruction, shows how it is associated with the hazard rate, and clarifies why it is life-threatening.

Directions for Future Research can Include
• Problem-specific derivations and specifications of the hazard-rate functions and ISODEs for the HRDI variables in various applications; • Implementation and improvement of the computational scenario formulated at any fixed s, t, and n x ∈  , is the conditional probability density of a random variable ( ) , and is such that Dynamic-Equilibrium Processes to the Survival Analyses Based on the Stochastic Quadratic-Hazard-Rate Model The remaining of the present work discusses a new application of the notion of DE DSP.In statistics, there are different quantitative approaches to survival analysis (e.g.[7]).Recently, interest has been in joint modeling, i.e. modeling longitudinal independent variables, or covariates, on top of time to event, whereas the covariates are modeled using a mixed-effects-model approach.The main idea of joint modeling is to complement conditional probability distributions, which are common in evolutions of multi-component populations and are conditioned with fixed values of deterministic parameters of the populations, with probability distributions of the component-individual parameters, thereby generalizing the deterministic parameters to their stochastic versions and taking into account their stochastic variability.Parameters that are not component-individual, remain deterministic.Joint modeling attracts growing attention in statistics.However, its main idea was developed in mathematical physics rather long ago (e.g.[8] [9]).It was motivated by needs in analysis of multi-component populations with a large number of components.This development treats the multi-componentness of the population as multi-modality of the probability density of stochastic parameters.It describes the parameters with the McKean-ISODE ([8], Section 6.2), which is even more general than ISODE (4) (not to mention linear Equation (18) below).
[3]the transition density ρ is stationary, i.e. depends on t s − only, rather than on both s and t, all processes in M ρ are homogeneous (e.g.[3], (2.2.8)).By the definition of the set M ρ , the marginal probability density of the process χ at time t is determined by the marginal density of the random of such densities shows that there can be a continuum of invariant processes that correspond to the same transition probability density.In other words, uniqueness is, in general, not fulfilled.More specifically, M ρ can contain a continuum of invariant processes, each of which corresponds to its individual invariant probability density.
If in M ρ there exists the DE process, it is of special importance in connection with convergence of processes that belong to M ρ .Convergence of One can show, under rather mild conditions, that if the above DE process exists, then processes that belong to M ρ converge in distribution, i.e. the marginal probability density of any of them at time t converges to the DE density Example 1, coincides with the invariant density.Thus, time is the DE process.Notably, since the expectation corresponding to DE density ( ) t y δ − , is t, the expectation of a DE process does not need be uniformly bounded in t.Example 4. Due to the well-known result (e.g.[3], (9.4.8) and limit relation (16), stationary density (14) is a DE density.It corresponds to the DSPs considered in Example 2. C t s is the Cauchy matrix of the ordinary differential

)
Expression (33) explicitly shows the following.Firstly, the bronchiole-individual hazard is the occlusive-layer thickness.Secondly, this hazard is obstructive because it, in the course of the age, grows.And, thirdly, this growth is life-threatening because the thickness, in the infinite-age limit, tends to the radius [27]eneral, scoring/grading systems should also take into account weight loss and muscle weakness, as well as the presence of other diseases.In agreement with this, a use of a number of non-negative and MF variables (including the GOLD grading) is in the core of clinical studies reported in[27]and other works.Also note that Remark 5 draws attention to the question on how specifically each entry of the HRDI-variable vector v can affect the bronchiole-individual hazard rate in its role indicated in (33).
[2]rnal of Applied Mathematics and Physics emphasized in Remark 4.However, in the problems of interest, the number of variables (or scalar equations) in system (4), n, is of the order of a few units or tens, or greater.How can one efficiently obtain probability distributions of such high dimensions from nonlinear ISODE (4)?Work[2]focuses on high-dimensional DSPs and ISODEs, and includes a few approaches to the problem.The simplest one provides an approximate model that comprises: As is above in Remark 6, Equation (4) is generally nonlinear in y.It is desirable that Equation (4) inherits the important advantage of linear Equation (18)