DENSERKS: Fortran sensitivity solvers using continuous, explicit Runge-Kutta schemes

DENSERKS is a Fortran sensitivity equation solver package designed for integrating models whose evolution can be described by ordinary differential equations (ODEs). A salient feature of DENSERKS is its support for both forward and adjoint sensitivity analyses, with built-in integrators for both first and second order continuous adjoint models. The software implements explicit Runge-Kutta methods with adaptive timestepping and high-order dense output schemes for the forward and the tangent linear model trajectory interpolation. Implementations of six Runge-Kutta methods are provided, with orders of accuracy ranging from two to eight. This makes DENSERKS suitable for a wide range of practical applications. The use of dense output, a novel approach in adjoint sensitivity analysis solvers, allows for a high-order cost-effective interpolation. This is a necessary feature when solving adjoints of nonlinear systems using highly accurate Runge-Kutta methods (order five and above). To minimize memory requirements and make long-time integrations computationally efficient, DENSERKS implements a two-level checkpointing mechanism. The code is tested on a selection of problems illustrating first and second order sensitivity analysis with respect to initial model conditions. The resulting derivative information is also used in a gradient-based optimization algorithm to minimize cost functionals dependent on a given set of model parameters

Discrete Second Order Adjoints in Atmospheric Chemical Transport Modeling

Atmospheric chemical transport models (CTMs) are essential tools for the study of air pollution, for environmental policy decisions, for the interpretation of observational data, and for producing air quality forecasts. Many air quality studies require sensitivity analyses, i.e., the computation of derivatives of the model output with respect to model parameters. The derivatives of a cost functional (defined on the model output) with respect to a large number of model parameters can be calculated efficiently through adjoint sensitivity analysis. While the traditional (first order) adjoint models give the gradient of the cost functional with respect to parameters, second order adjoint models give second derivative information in the form of products between the Hessian of the cost functional and a user defined vector. In this paper we discuss the mathematical foundations of the discrete second order adjoint sensitivity method and present a complete set of computational tools for performing second order sensitivity studies in three-dimensional atmospheric CTMs. The tools include discrete second order adjoints of Runge Kutta and of Rosenbrock time stepping methods for stiff equations together with efficient implementation strategies. Numerical examples illustrate the use of these computational tools in important applications like sensitivity analysis, optimization, uncertainty quantification, and the calculation of directions of maximal error growth in three-dimensional atmospheric CTMs

Experiment design for systems biology

Thesis (Ph. D.)--Massachusetts Institute of Technology, Dept. of Biological Engineering, 2009.Cataloged from PDF version of thesis.Includes bibliographical references (p. 219-233).Mechanism-based chemical kinetic models are increasingly being used to describe biological signaling. Such models serve to encapsulate current understanding of pathways and to enable insight into complex biological processes. Despite the growing interest in these models, a number of challenges frustrate the construction of high-quality models. First, the chemical reactions that control biochemical processes are only partially known, and multiple, mechanistically distinct models often fit all of the available data and known chemistry. We address this by providing methods for designing dynamic stimuli that can distinguish among models with different reaction mechanisms in stimulus-response experiments. We evaluated our method on models of antibody-ligand binding, mitogen-activated protein kinase phosphorylation and de-phosphorylation, and larger models of the epidermal growth factor receptor (EGFR) pathway. Inspired by these computational results, we tested the idea that pulses of EGF could help elucidate the relative contribution of different feedback loops within the EGFR network. These experimental results suggest that models from the literature do not accurately represent the relative strength of the various feedback loops in this pathway. In particular, we observed that the endocytosis and feedback loop was less strong than predicted by models, and that other feedback mechanisms were likely necessary to deactivate ERK after EGF stimulation. Second, chemical kinetic models contain many unknown parameters, at least some of which must be estimated by fitting to time-course data. We examined this question in the context of a pathway model of EGF and neuronal growth factor (NGF) signaling. Computationally, we generated a palette of experimental perturbation data that included different doses of EGF and NGF as well as single and multiple gene knockdowns and overexpressions. While no single experiment could accurately estimate all of the parameters, we identified a set of five complementary experiments that could. These results suggest that there is reason to be optimistic about the prospects for parameter estimation in even large models. Third, there is no standard formulation for chemical kinetic models of biological signaling. We propose a general and concise formulation of mass action kinetics based on sparse matrices and Kronecker products. This formulation allows any mass action model and its partial derivatives to be represented by simple matrix equations, which enabled straightforward application of several numerical methods. We show that models that use other rate laws such as MichaelisMenten can be converted to our formulation. We demonstrate this by converting a model of Escherichia coli central carbon metabolism to use only mass action kinetics. The dynamics of the new model are similar to the original model. However, we argue that because our model is based on fewer approximations it has the potential to be more accurate over a wider range of conditions. Taken together, the work presented here demonstrates that experimental design methodology can be successfully used to improve the quality of mechanism-based chemical kinetic models.by Joshua Farley Apgar.Ph.D

Discrete Adjoints: Theoretical Analysis, Efficient Computation, and Applications

The technique of automatic differentiation provides directional derivatives and discrete adjoints with working accuracy. A complete complexity analysis of the basic modes of automatic differentiation is available. Therefore, the research activities are focused now on different aspects of the derivative calculation, as for example the efficient implementation by exploitation of structural information, studies of the theoretical properties of the provided derivatives in the context of optimization problems, and the development and analysis of new mathematical algorithms based on discrete adjoint information. According to this motivation, this habilitation presents an analysis of different checkpointing strategies to reduce the memory requirement of the discrete adjoint computation. Additionally, a new algorithm for computing sparse Hessian matrices is presented including a complexity analysis and a report on practical experiments. Hence, the first two contributions of this thesis are dedicated to an efficient computation of discrete adjoints. The analysis of discrete adjoints with respect to their theoretical properties is another important research topic. The third and fourth contribution of this thesis focus on the relation of discrete adjoint information and continuous adjoint information for optimal control problems. Here, differences resulting from different discretization strategies as well as convergence properties of the discrete adjoints are analyzed comprehensively. In the fifth contribution, checkpointing approaches that are successfully applied for the computation of discrete adjoints, are adapted such that they can be used also for the computation of continuous adjoints. Additionally, the fifth contributions presents a new proof of optimality for the binomial checkpointing that is based on new theoretical results. Discrete adjoint information can be applied for example for the approximation of dense Jacobian matrices. The development and analysis of new mathematical algorithms based on these approximate Jacobians is the topic of the sixth contribution. Is was possible to show global convergence to first-order critical points for a whole class of trust-region methods. Here, the usage of inexact Jacobian matrices allows a considerable reduction of the computational complexity

Sensitivity analysis of oscillating dynamical systems with applications to the mammalian circadian clock

Thesis (Ph. D.)--Massachusetts Institute of Technology, Dept. of Chemical Engineering, 2008.This electronic version was submitted by the student author. The certified thesis is available in the Institute Archives and Special Collections.Includes bibliographical references (p. 227-234).The work presented in this thesis consists of two major parts. In Chapter 2, the theory for sensitivity analysis of oscillatory systems is developed and discussed. Several contributions are made, in particular in the precise definition of phase sensitivities and in the generalization of the theory to all types of autonomous oscillators. All methods rely on the solution of a boundary value problem, which identifies the periodic orbit. The choice of initial condition on the limit cycle has important consequences for phase sensitivity analysis, and its influence is quantified and discussed in detail. The results are exact and efficient to compute compared to existing partial methods. The theory is then applied to different models of the mammalian circadian clock system in the following chapters. First, different types of sensitivities in a pair of smaller models are analyzed. The models have slightly different architectures, with one having an additional negative feedback loop compared to the other. The differences in their behavior with respect to phases, the period and amplitude are discussed in the context of their network architecture. It is found that, contrary to previous assumptions in the literature, the additional negative feedback loop makes the model less "flexible" in at least one sense that was studied here. The theory was also applied to larger, more detailed models of the mammalian circadian clock, based on the original model of Forger and Peskin. Between the original model's publication in 2003 and the present time, several key advances were made in understanding the mechanistic detail of the mammalian circadian clock, and at least one additional clock gene was identified. These advances are incorporated in an extended model, which is then studied using sensitivity analysis. Period sensitivity analysis is performed first and it was found that only one negative feedback loop dominates the setting of the period.(cont.) This was an interesting one-to-one correlation between one topological feature of the network and a single metric of network performance. This led to the question of whether the network architecture is modular, in the sense that each of the several feedback loops might be responsible for a separate network function. A function of particular interest is the ability to separately track "dawn" and "dusk", which is reported to be present in the circadian clock. The ability of the mammalian circadian clock to modify different relative phases --defined by different molecular events -- independently of the period was analyzed. If the model can maintain a perceived day -- defined by the time difference between two phases -- of different lengths, it can be argued that the model can track dawn and dusk separately. This capability is found in all mammalian clock models that were studied in this work, and furthermore, that a network-wide effort is needed to do so. Unlike in the case of the period sensitivities, relative phase sensitivities are distributed throughout several feedback loops. Interestingly, a small number of "key parameters" could be identified in the detailed models that consistently play important roles in the setting of period, amplitude and phases. It appears that most circadian clock features are under shared control by local parameters and by the more global "key parameters". Lastly, it is shown that sensitivity analysis, in particular period sensitivity analysis, can be very useful in parameter estimation for oscillatory systems biology models. In an approach termed "feature-based parameter fitting", the model's parameter values are selected based on their impact on the "features" of an oscillation (period, phases, amplitudes) rather than concentration data points. It is discussed how this approach changes the cost function during the parameter estimation optimization, and when it can be beneficial.(cont.) A minimal model system from circadian biology, the Goodwin oscillator, is taken as an example. Overall, in this thesis it is shown that the contributions made to the theoretical understanding of sensitivities in oscillatory systems are relevant and useful in trying to answer questions that are currently open in circadian biology. In some cases, the theory could indicate exactly which experiments or detailed mechanistic studies are needed in order to perform meaningful mathematical analysis of the system as a whole. It is shown that, provided the biologically relevant quantities are analyzed, a network-wide understanding of the interplay between network function and topology can be gained and differences in performance between models of different size or topology can be quantified.by Anna Katharina Wilkins.Ph.D

Nonsmooth dynamic optimization of systems with varying structure

Thesis (Ph. D.)--Massachusetts Institute of Technology, Dept. of Mechanical Engineering, 2011.Cataloged from PDF version of thesis.Includes bibliographical references (p. 357-365).In this thesis, an open-loop numerical dynamic optimization method for a class of dynamic systems is developed. The structure of the governing equations of the systems under consideration change depending on the values of the states, parameters and the controls. Therefore, these systems are called systems with varying structure. Such systems occur frequently in the models of electric and hydraulic circuits, chemical processes, biological networks and machinery. As a result, the determination of parameters and controls resulting in the optimal performance of these systems has been an important research topic. Unlike dynamic optimization problems where the structure of the underlying system is constant, the dynamic optimization of systems with varying structure requires the determination of the optimal evolution of the system structure in time in addition to optimal parameters and controls. The underlying varying structure results in nonsmooth and discontinuous optimization problems. The nonsmooth single shooting method introduced in this thesis uses concepts from nonsmooth analysis and nonsmooth optimization to solve dynamic optimization problems involving systems with varying structure whose dynamics can be described by locally Lipschitz continuous ordinary or differential-algebraic equations. The method converts the infinitedimensional dynamic optimization problem into an nonlinear program by parameterizing the controls. Unlike the state of the art, the method does not enumerate possible structures explicitly in the optimization and it does not depend on the discretization of the dynamics. Instead, it uses a special integration algorithm to compute state trajectories and derivative information. As a result, the method produces more accurate solutions to problems where the underlying dynamics is highly nonlinear and/or stiff for less effort than the state of the art. The thesis develops substitutes for the gradient and the Jacobian of a function in case these quantities do not exist. These substitutes are set-valued maps and an elements of these maps need to be computed for optimization purposes. Differential equations are derived whose solutions furnish the necessary elements. These differential equations have discontinuities in time. A numerical method for their solution is proposed based on state event location algorithms that detects these discontinuities. Necessary conditions of optimality for nonlinear programs are derived using these substitutes and it is shown that nonsmooth optimization methods called bundle methods can be used to obtain solutions satisfying these necessary conditions. Case studies compare the method to the state of the art and investigate its complexity empirically.by Mehmet Yunt.Ph.D