Stochastic Reaction Networks Within Interacting Compartments.

Stochastic reaction networks, which are usually modeled as continuous-time Markov chains on [Formula: see text], and simulated via a version of the "Gillespie algorithm," have proven to be a useful tool for the understanding of processes, chemical and otherwise, in homogeneous environments. There are multiple avenues for generalizing away from the assumption that the environment is homogeneous, with the proper modeling choice dependent upon the context of the problem being considered. One such generalization was recently introduced in Duso and Zechner (Proc Nat Acad Sci 117(37):22674-22683 , Duso and Zechner (2020)), where the proposed model includes a varying number of interacting compartments, or cells, each of which contains an evolving copy of the stochastic reaction system. The novelty of the model is that these compartments also interact via the merging of two compartments (including their contents), the splitting of one compartment into two, and the appearance and destruction of compartments. In this paper we begin a systematic exploration of the mathematical properties of this model. We (i) obtain basic/foundational results pertaining to explosivity, transience, recurrence, and positive recurrence of the model, (ii) explore a number of examples demonstrating some possible non-intuitive behaviors of the model, and (iii) identify the limiting distribution of the model in a special case that generalizes three formulas from an example in Duso and Zechner (Proc Nat Acad Sci 117(37):22674-22683 , Duso and Zechner (2020)).

Time-dependent product-form Poisson distributions for reaction networks with higher order complexes.

It is well known that stochastically modeled reaction networks that are complex balanced admit a stationary distribution that is a product of Poisson distributions. In this paper, we consider the following related question: supposing that the initial distribution of a stochastically modeled reaction network is a product of Poissons, under what conditions will the distribution remain a product of Poissons for all time? By drawing inspiration from Crispin Gardiner's "Poisson representation" for the solution to the chemical master equation, we provide a necessary and sufficient condition for such a product-form distribution to hold for all time. Interestingly, the condition is a dynamical "complex-balancing" for only those complexes that have multiplicity greater than or equal to two (i.e. the higher order complexes that yield non-linear terms to the dynamics). We term this new condition the "dynamical and restricted complex balance" condition (DR for short).

Low Variance Couplings for Stochastic Models of Intracellular Processes with Time-Dependent Rate Functions.

A number of coupling strategies are presented for stochastically modeled biochemical processes with time-dependent parameters. In particular, the stacked coupling is introduced and is shown via a number of examples to provide an exceptionally low variance between the generated paths. This coupling will be useful in the numerical computation of parametric sensitivities and the fast estimation of expectations via multilevel Monte Carlo methods. We provide the requisite estimators in both cases.

Discrepancies between extinction events and boundary equilibria in reaction networks.

Reaction networks are mathematical models of interacting chemical species that are primarily used in biochemistry. There are two modeling regimes that are typically used, one of which is deterministic and one that is stochastic. In particular, the deterministic model consists of an autonomous system of differential equations, whereas the stochastic system is a continuous-time Markov chain. Connections between the two modeling regimes have been studied since the seminal paper by Kurtz (J Chem Phys 57(7):2976-2978, 1972), where the deterministic model is shown to be a limit of a properly rescaled stochastic model over compact time intervals. Further, more recent studies have connected the long-term behaviors of the two models when the reaction network satisfies certain graphical properties, such as weak reversibility and a deficiency of zero. These connections have led some to conjecture a link between the long-term behavior of the two models exists, in some sense. In particular, one is tempted to believe that positive recurrence of all states for the stochastic model implies the existence of positive equilibria in the deterministic setting, and that boundary equilibria of the deterministic model imply the occurrence of an extinction event in the stochastic setting. We prove in this paper that these implications do not hold in general, even if restricting the analysis to networks that are bimolecular and that conserve the total mass. In particular, we disprove the implications in the special case of models that have absolute concentration robustness, thus answering in the negative a conjecture stated in the literature in 2014.

Non-explosivity of Stochastically Modeled Reaction Networks that are Complex Balanced.

We consider stochastically modeled reaction networks and prove that if a constant solution to the Kolmogorov forward equation decays fast enough relatively to the transition rates, then the model is non-explosive. In particular, complex-balanced reaction networks are non-explosive.

Conditions for extinction events in chemical reaction networks with discrete state spaces.

We study chemical reaction networks with discrete state spaces and present sufficient conditions on the structure of the network that guarantee the system exhibits an extinction event. The conditions we derive involve creating a modified chemical reaction network called a domination-expanded reaction network and then checking properties of this network. Unlike previous results, our analysis allows algorithmic implementation via systems of equalities and inequalities and suggests sequences of reactions which may lead to extinction events. We apply the results to several networks including an EnvZ-OmpR signaling pathway in Escherichia coli.

Product-Form Stationary Distributions for Deficiency Zero Networks with Non-mass Action Kinetics.

In many applications, for example when computing statistics of fast subsystems in a multiscale setting, we wish to find the stationary distributions of systems of continuous-time Markov chains. Here we present a class of models that appears naturally in certain averaging approaches whose stationary distributions can be computed explicitly. In particular, we study continuous-time Markov chain models for biochemical interaction systems with non-mass action kinetics whose network satisfies a certain constraint. Analogous with previous related results, the distributions can be written in product form.

Stochastic representations of ion channel kinetics and exact stochastic simulation of neuronal dynamics.

In this paper we provide two representations for stochastic ion channel kinetics, and compare the performance of exact simulation with a commonly used numerical approximation strategy. The first representation we present is a random time change representation, popularized by Thomas Kurtz, with the second being analogous to a "Gillespie" representation. Exact stochastic algorithms are provided for the different representations, which are preferable to either (a) fixed time step or (b) piecewise constant propensity algorithms, which still appear in the literature. As examples, we provide versions of the exact algorithms for the Morris-Lecar conductance based model, and detail the error induced, both in a weak and a strong sense, by the use of approximate algorithms on this model. We include ready-to-use implementations of the random time change algorithm in both XPP and Matlab. Finally, through the consideration of parametric sensitivity analysis, we show how the representations presented here are useful in the development of further computational methods. The general representations and simulation strategies provided here are known in other parts of the sciences, but less so in the present setting.

Lyapunov Functions, Stationary Distributions, and Non-equilibrium Potential for Reaction Networks.

We consider the relationship between stationary distributions for stochastic models of reaction systems and Lyapunov functions for their deterministic counterparts. Specifically, we derive the well-known Lyapunov function of reaction network theory as a scaling limit of the non-equilibrium potential of the stationary distribution of stochastically modeled complex balanced systems. We extend this result to general birth-death models and demonstrate via example that similar scaling limits can yield Lyapunov functions even for models that are not complex or detailed balanced, and may even have multiple equilibria.

Hybrid pathwise sensitivity methods for discrete stochastic models of chemical reaction systems.

Stochastic models are often used to help understand the behavior of intracellular biochemical processes. The most common such models are continuous time Markov chains (CTMCs). Parametric sensitivities, which are derivatives of expectations of model output quantities with respect to model parameters, are useful in this setting for a variety of applications. In this paper, we introduce a class of hybrid pathwise differentiation methods for the numerical estimation of parametric sensitivities. The new hybrid methods combine elements from the three main classes of procedures for sensitivity estimation and have a number of desirable qualities. First, the new methods are unbiased for a broad class of problems. Second, the methods are applicable to nearly any physically relevant biochemical CTMC model. Third, and as we demonstrate on several numerical examples, the new methods are quite efficient, particularly if one wishes to estimate the full gradient of parametric sensitivities. The methods are rather intuitive and utilize the multilevel Monte Carlo philosophy of splitting an expectation into separate parts and handling each in an efficient manner.

Comparison of finite difference based methods to obtain sensitivities of stochastic chemical kinetic models.

Sensitivity analysis is a powerful tool in determining parameters to which the system output is most responsive, in assessing robustness of the system to extreme circumstances or unusual environmental conditions, in identifying rate limiting pathways as a candidate for drug delivery, and in parameter estimation for calculating the Hessian of the objective function. Anderson [SIAM J. Numer. Anal. 50, 2237 (2012)] shows the advantages of the newly developed coupled finite difference (CFD) estimator over the common reaction path (CRP) [M. Rathinam, P. W. Sheppard, and M. Khammash, J. Chem. Phys. 132, 034103 (2010)] estimator. In this paper, we demonstrate the superiority of the CFD estimator over the common random number (CRN) estimator in a number of scenarios not considered previously in the literature, including the sensitivity of a negative log likelihood function for parameter estimation, the sensitivity of being in a rare state, and a sensitivity with fast fluctuating species. In all examples considered, the superiority of CFD over CRN is demonstrated. We also provide an example in which the CRN method is superior to the CRP method, something not previously observed in the literature. These examples, along with Anderson's results, lead to the conclusion that CFD is currently the best estimator in the class of finite difference estimators of stochastic chemical kinetic models.

Mixing times for two classes of stochastically modeled reaction networks.

The past few decades have seen robust research on questions regarding the existence, form, and properties of stationary distributions of stochastically modeled reaction networks. When a stochastic model admits a stationary distribution an important practical question is: what is the rate of convergence of the distribution of the process to the stationary distribution? With the exception of [1] pertaining to models whose state space is restricted to the non-negative integers, there has been a notable lack of results related to this rate of convergence in the reaction network literature. This paper begins the process of filling that hole in our understanding. In this paper, we characterize this rate of convergence, via the mixing times of the processes, for two classes of stochastically modeled reaction networks. Specifically, by applying a Foster-Lyapunov criteria we establish exponential ergodicity for two classes of reaction networks introduced in [2]. Moreover, we show that for one of the classes the convergence is uniform over the initial state.

Simultaneous tectonic Descemet stripping endothelial keratoplasty and tectonic Bowman layer transplant for the management of corneal perforation.

AIM: To report the outcomes of three cases of corneal perforation managed with simultaneous tectonic Descemet stripping endothelial keratoplasty (t-DSEK) and tectonic Bowman layer transplant (t-BLT) as an alternative to tectonic penetrating keratoplasty (t-PKP). METHODS: Three eyes of three patients receiving simultaneous t-DSEK and t-BLT for corneal perforation were included. The technique for DSEK was modified depending on individual requirements. The t-BLT technique was standardised using an 8 mm graft and fixated with a running suture. Success was measured by the ability of this procedure to close a corneal perforation. RESULTS: All three cases achieved tectonic eye globe restoration and remained stable during the minimum 3-month observation period. Reinterventions were relatively common: 2 cases required amniotic membrane transplant for persistent epithelial defects. One case required DSEK rebubbling. One case developed angle closure glaucoma requiring surgical peripheral iridectomy. CONCLUSION: Simultaneous t-DSEK and t-BLT may be a useful strategy for the management of corneal perforation as an alternative management to t-PKP for selected cases.

A finite difference method for estimating second order parameter sensitivities of discrete stochastic chemical reaction networks.

We present an efficient finite difference method for the approximation of second derivatives, with respect to system parameters, of expectations for a class of discrete stochastic chemical reaction networks. The method uses a coupling of the perturbed processes that yields a much lower variance than existing methods, thereby drastically lowering the computational complexity required to solve a given problem. Further, the method is simple to implement and will also prove useful in any setting in which continuous time Markov chains are used to model dynamics, such as population processes. We expect the new method to be useful in the context of optimization algorithms that require knowledge of the Hessian.

Case series of staged lamellar corneal grafting with thin manual Descemet's stripping endothelial keratoplasty, followed by manual deep anterior lamellar keratoplasty, as an alternative to penetrating keratoplasty.

The purpose of this report is to present the outcomes of Descemet's stripping endothelial keratoplasty (DSEK), followed by deep anterior lamellar keratoplasty (DALK), as an alternative to penetrating keratoplasty (PKP) for different indications. Patients in this retrospective case series underwent manual DSEK, followed by manual DALK, for tectonic and/or visual reasons. It includes three cases that underwent DSEK followed by DALK as an alternative to PKP for the following conditions: failed PKP and stromal scarring, aphakic bullous keratopathy and herpetic stromal scarring, and herpetic corneal perforation. All cases had successful anatomical repair with visual improvement. Many corneal conditions involving endothelium and stroma have been traditionally managed with PKP. However, in certain settings with high risk of complications, this might not be an appropriate option. Therefore, for these cases, we propose management with DSEK followed by DALK as an alternative approach in a 'closed anterior chamber' fashion, which can reduce the risk of significant, sometimes irreversible, complications.

One-Year Front Versus Central and Paracentral Corneal Changes After Bowman Layer Transplantation for Keratoconus.

PURPOSE: To report the front corneal versus central and paracentral corneal changes after Bowman layer transplantation for keratoconus in a tertiary hospital in the United Kingdom. METHODS: Five eyes of 5 patients receiving Bowman layer transplant for advanced keratoconus in Royal Gwent Hospital (Newport, United Kingdom) were included. Preoperative and postoperative visual acuity; Kmax; Kmean, and corneal cylinder in the front cornea, 4.5 mm central, and 6 mm central; and corneal thickness were analyzed. RESULTS: Corneal flattening and reduction in corneal astigmatism was observed, more marked in the central and paracentral zone, allowing for improvement in best-corrected visual acuity with the aid of visual correction in 4 eyes. CONCLUSIONS: These results support previous data reporting Bowman layer transplantation as a useful strategy in the treatment of advanced keratoconus and suggest greater attention may be focused on central or paracentral corneal changes.

1 year posterior corneal changes after Bowman Layer Transplant for keratoconus.

INTRODUCTION: To report the posterior corneal changes after Bowman Layer Transplant for keratoconus in a tertiary hospital in the UK. METHODS: 5 eyes of 5 patients receiving Bowman Layer Transplant for advanced keratoconus in Royal Gwent Hospital (Newport, UK) were included. Pre and postoperative posterior corneal astigmatism, posterior Kmean, and back surface elevation were analysed. RESULTS: No significant changes were seen in the posterior corneal astigmatism, posterior Kmean, or back surface elevation between the pre- and postoperative period. CONCLUSION: This results would support the idea that the corneal changes seen after Bowman Layer Transplant are mainly in the anterior corneal surface.

On reaction network implementations of neural networks.

This paper is concerned with the utilization of deterministically modelled chemical reaction networks for the implementation of (feed-forward) neural networks. We develop a general mathematical framework and prove that the ordinary differential equations (ODEs) associated with certain reaction network implementations of neural networks have desirable properties including (i) existence of unique positive fixed points that are smooth in the parameters of the model (necessary for gradient descent) and (ii) fast convergence to the fixed point regardless of initial condition (necessary for efficient implementation). We do so by first making a connection between neural networks and fixed points for systems of ODEs, and then by constructing reaction networks with the correct associated set of ODEs. We demonstrate the theory by constructing a reaction network that implements a neural network with a smoothed ReLU activation function, though we also demonstrate how to generalize the construction to allow for other activation functions (each with the desirable properties listed previously). As there are multiple types of 'networks' used in this paper, we also give a careful introduction to both reaction networks and neural networks, in order to disambiguate the overlapping vocabulary in the two settings and to clearly highlight the role of each network's properties.

Product-form stationary distributions for deficiency zero chemical reaction networks.

We consider stochastically modeled chemical reaction systems with mass-action kinetics and prove that a product-form stationary distribution exists for each closed, irreducible subset of the state space if an analogous deterministically modeled system with mass-action kinetics admits a complex balanced equilibrium. Feinberg's deficiency zero theorem then implies that such a distribution exists so long as the corresponding chemical network is weakly reversible and has a deficiency of zero. The main parameter of the stationary distribution for the stochastically modeled system is a complex balanced equilibrium value for the corresponding deterministically modeled system. We also generalize our main result to some non-mass-action kinetics.

Tectonic Descemet Stripping Endothelial Keratoplasty for the Management of Corneal Perforation: A Case Series.

PURPOSE: Corneal perforations pose a considerable challenge for ophthalmic surgeons. The aim remains restoring the anatomical integrity of eye while attempting to preserve as much visual function as possible. To our knowledge, we report the first case series of 4 successful tectonic Descemet stripping endothelial keratoplasty (DSEK) grafts completed for acute corneal perforation. In all cases, restoration of globe integrity was achieved. In cases where visual potential remained, deep anterior lamellar keratoplasty over DSEK was offered. METHODS: Four patients presenting to Royal Gwent Hospital (Newport, Wales) with corneal perforation were included. Etiological origins included exposure keratopathy, corneal hydrops secondary to pellucid marginal degeneration with fistulization-associated aqueous leak, and herpetic keratitis. Patients were treated with initial temporary therapies, including glue patch, the use of bandage contact lenses, amniotic membrane transplant, and antibiotic cover, depending on the size and location of the perforation. Subsequent DSEK was subsequently performed under local anesthesia. RESULTS: All patients had successful tectonic grafts (3 DSEKs and 1 hemi-DSEK) up to 4 months postoperatively. In the case involving a hemi-DSEK, the patient underwent 2 refloating attempts of the graft before anchoring suture placement because of repeated graft detachment. Successful graft placement and tectonic globe restoration was confirmed by anterior segment imaging. CONCLUSIONS: Tectonic DSEK provides corneal surgeons with a new modality of treating corneal perforations. It provides clinicians an additional tool in their armamentarium in complex cases where anterior lamellar or penetrating keratoplasty may be contraindicated or deemed high risk of complications.