Search

There is a continuous demand for new and improved methods of understanding our climate system. The work in this thesis focuses on the study of delayed feedback and critical transitions. There is much room to develop upon these concepts in their application to the climate system. We explore the two concepts independently, but also note that the two are not mutually exclusive. The thesis begins with a review of delay differential equation (DDE) theory and the use of delay models in climate, followed by a review of the literature on critical transitions and examples of critical transitions in climate. We introduce various methods of deriving delay models from more complex systems. Our main results center around the Saltzman and Maasch (1988) model for the Pleistocene climate (`Carbon cycle instability as a cause of the late Pleistocene ice age oscillations: modelling the asymmetric response.' Global biogeochemical cycles, 2(2):177-185, 1988). We observe that the model contains a chain of first-order reactions. Feedback chains of this type limits to a discrete delay for long chains. We can then approximate the chain by a delay, resulting in scalar DDE for ice mass. Through bifurcation analysis under varying the delay, we discover a previously unexplored bistable region and consider solutions in this parameter region when subjected to periodic and astronomical forcing. The astronomical forcing is highly quasiperiodic, containing many overlapping frequencies from variations in the Earth's orbit. We find that under the astronomical forcing, the model exhibits a transition in time that resembles what is seen in paleoclimate records, known as the Mid-Pleistocene Transition. This transition is a distinct feature of the quasiperiodic forcing, as confi rmed by the change in sign of the leading nite-time Lyapunov exponent. Additional results involve a box model of the Atlantic meridional overturning circulation under a future climate scenario and time-dependent freshwater forcing. We find that the model exhibits multiple types of critical transitions, as well as recovery from potential critical transitions. We conclude with an outlook on how the work presented in this thesis can be utilised for further studies of the climate system and beyond.

This thesis is concerned with mathematical modelling of oncolytic virotherapy: the use of genetically modified viruses to selectively spread, replicate and destroy cancerous cells in solid tumours. Traditional spatially-dependent modelling approaches have previously assumed that virus spread is due to viral diffusion in solid tumours, and also neglect the time delay introduced by the lytic cycle for viral replication within host cells. A deterministic, age-structured reaction-diffusion model is developed for the spatially-dependent interactions of uninfected cells, infected cells and virus particles, with the spread of virus particles facilitated by infected cell motility and delay. Evidence of travelling wave behaviour is shown, and an asymptotic approximation for the wave speed is derived as a function of key parameters. Next, the same physical assumptions as in the continuum model are used to develop an equivalent discrete, probabilistic model for that is valid in the limit of low particle concentrations. This mesoscopic, compartment-based model is then validated against known test cases, and it is shown that the localised nature of infected cell bursts leads to inconsistencies between the discrete and continuum models. The qualitative behaviour of this stochastic model is then analysed for a range of key experimentally-controllable parameters. Two-dimensional simulations of in vivo and in vitro therapies are then analysed to determine the effects of virus burst size, length of lytic cycle, infected cell motility, and initial viral distribution on the wave speed, consistency of results and overall success of therapy. Finally, the experimental difficulty of measuring the effective motility of cells is addressed by considering effective medium approximations of diffusion through heterogeneous tumours. Considering an idealised tumour consisting of periodic obstacles in free space, a two-scale homogenisation technique is used to show the effects of obstacle shape on the effective diffusivity. A novel method for calculating the effective continuum behaviour of random walks on lattices is then developed for the limiting case where microscopic interactions are discrete.

This thesis concerns the development of a method for the detection of small solutions to delay differential equations. The detection of small solutions is important because their presence has significant influence on the analytical prop¬erties of an equation. However, to date, analytical methods are of only limited practical use. Therefore this thesis focuses on the development of a reliable new method, based on finite order approximations of the underlying infinite dimen¬sional problem, which can detect small solutions. Decisions (concerning the existence, or otherwise, of small solutions) based on our visualisation technique require an understanding of the underlying methodol¬ogy behind our approach. Removing this need would be attractive. The method we have developed can be automated, and at the end of the thesis we present a prototype Matlab code for the automatic detection of small solutions to delay differential equations.

In many service systems, customers are often presented with information about queue lengths or waiting times. This information has the potential to influence the decisions of customers on whether or not to join a queue to receive service. However, it is often the case that the information that customers receive is not given in real time and has some inherent delay. For example, the information provided could be updated periodically and therefore the customer receives information about the system from some time in the past. In this thesis, we focus our attention on fluid models of queueing systems. Generally, the fluid models in this thesis are the limiting objects of scaled stochastic systems and take the form of delay differential equations. We show that the dynamics of the queueing system can change dramatically depending on how large the delay in information is. Thus, our main goal in this thesis is to explore and understand how the type of delay and the size of the delay impact the underlying dynamics of different queueing models.

Many service systems use internet or smartphone app technology to notify customers about their expected waiting times or queue lengths via delay announcements, allowing customers to decide which queue to join. However, in many cases, either the information might be delayed or customers might require time to travel to the queue of their choice, thus causing a lag in information. We model multiple-queue systems through delay differential equations, and study how the delay in information affects the dynamics of the queues. When the delay is sufficiently large, the queues may oscillate indefinitely throughout time. We develop accurate approximations for the amplitude of these oscillations by implementing two numerical methods. The first technique is a classical analytic method that yields a closed-form approximation in terms of the model parameters. The second approximation method is novel, and it uses a statistical technique to deliver highly accurate approximations over a wider range of parameters. The oscillations in queue lengths are generally undesirable both for the service providers and the customers. This motivates us to explore how the delay announcement can be used to limit the oscillations. We show that, in some cases, using information about queue's velocity (the rate at which the queue length is changing) in the delay announcement can eliminate oscillations created by delays in information. We derive a fixed point equation for determining the optimal amount of velocity information that should be used and find closed form upper and lower bounds on its value. When the oscillations cannot be eliminated altogether, we identify the amount of velocity information that minimizes the amplitude of the oscillations. However, we also find that using too much velocity information can create oscillations in the queue lengths that would otherwise be stable. When the delay in information is caused by customers traveling to the queues, the delay may vary from customer to customer. We propose a queueing model that treats the individual's delay not as a constant, but as a random variable drawn from a fixed distribution. This generalized model allows us to identify the properties of dynamics that are independent of the delay distribution, such as the existence and uniqueness of the equilibrium state. However, the stability of the equilibrium and the presence of oscillations depend on the delay distribution. We therefore give an overview of the system's stability region for different common delay distributions, and finally offer a numerical method of approximating the stability region when the distribution is unknown.

This thesis considers the problem of robust geometric regulation for tracking and disturbance rejection of systems governed by delay differential equations. It is well known that geometric regulation can be highly sensitive to system parameters and hence such designs are not always robust. In particular, when employing numerical approximations to delay systems, the resulting finite dimensional models inherit natural approximation errors that can impact robustness. This demonstrates this lack of robustness and then addresses robustness by employing versions of robust regulation that have been developed for infinite dimensional systems. Numerical examples are given to illustrate the ideas and to test the robustness of the regulator.

Delay differential equations (DDEs) are often used to model systems with time-delayed effects, and they have found applications in fields such as climate dynamics, biosciences, engineering, and control theory. In contrast to ordinary differential equations (ODEs), the phase space associated even with a scalar DDE is infinite-dimensional. Oftentimes, it is desirable to have low-dimensional ODE systems that capture qualitative features as well as approximate certain quantitative aspects of the DDE dynamics. In this thesis, we present a Galerkin scheme for a broad class of DDEs and derive convergence results for this scheme. In contrast to other Galerkin schemes devised in the DDE literature, the main new ingredient here is the use of the so called Koornwinder polynomials, which are orthogonal polynomials under an inner product with a point mass. A main advantage of using such polynomials is that they live in the domain of the underlying linear operator, which arguably simplifies the related numerical treatments. The obtained results generalize a previous work to the case of DDEs with multiply delays in the linear terms, either discrete or distributed, or both. We also consider the more challenging case of discrete delays in the nonlinearity and obtain a convergence result by assuming additional assumptions about the Galerkin approximations of the linearized systems.

The standard eigenvalue problem of finding the zeros of det(zI − A) is ubiquitous, and comes from studying solutions to x′ = Ax as well as myriad other sources. Similarly, if x′(t) = Ax(t) + Bx(t − 1) is some delay-differential equation (arising, say, from modeling the spread of a disease, or from population growth models), then stability is determined by computing the roots of det(zI − A − Be−z). Models of physical processes where energy slowly leaks away to infinity lead to similar problems. These physical systems are typically modeled in terms of differential equations which are then discretized using e.g. collocation or finite elements. For example, such a discretization gives a correspondence between the quantum scattering resonances associated to a quantum corral and the zeros of det(A − zB + C(z)), where A, B, C(z) ∈ ℂn×n and the highly nonlinear entries of C(z) involve square roots and ratios of Bessel functions. In each case, we are led to so-called nonlinear eigenvalue problems of the form T(λ)v = 0, v ≠ 0, where T : Ω ⊂ ℂ → ℂn×n is a matrix-valued function, and λ is called an eigenvalue of T. The first contribution of this thesis is theorems for localizing eigenvalues of general matrix-valued functions, effectively reducing the region in which eigenvalues of T are known to lie from all of Ω down to a smaller space, and deducing eigenvalue counts within regions that meet certain conditions. These theorems are derived by working with the diagonal entries or diagonal blocks of T, such as our generalization of Gershgorin’s theorem, or by considering nonlinear generalizations of pseudospectra. Localization and counting results allow better initial guesses or shifts for iterative algorithms, guide the selection of an appropriate closed contour for contour integral-based algorithms, and facilitate error analysis in cases where eigenvalues can be confined to tiny regions. The second contribution of this thesis is to exploit these results in several contexts. To start with, a variety of strategies for getting the most out of our main localization theorems will be presented and applied to several test problems. Then we foray into the analysis of higher-order and delay differential equations, using our localization results to help bound asymptotic growth of solutions, and using our generalization of the notion of pseudospectra to concretely bound transient growth both above and below; a model for a semiconductor laser with phase-conjugate feedback acts as the central example. The last application we will treat is about the resonances for electrons trapped in circular quantum corrals, microscopic structures built by placing atoms in a circle on a metal surface. We provide a framework for comparing various elastic scattering models, and use it to bound the error between resonances computed from the naïve particle-in-a-box model and resonances computed from a model that takes quantum tunneling into account.

With the emergence of new measurement devices, the non-smooth nature of the borehole geometry has been comprehended more accurately. In many happenstances, the borehole has been seen taking the shape of a corkscrew or helix. Referred to as micro-tortuosity, or more commonly spiraling, this drilling dysfunction correlates with lower-than-expected rates of penetration, increased shocks and vibrations, damage to components and tools, and smaller drift diameters. In this dissertation, the relevant mechanism and parameters leading to borehole spiraling are identified and studied. Spiraling is predicted by conducting a stability analysis of the linearized delay differential equations governing the borehole propagation. These evolution equations, expressed in terms of the borehole inclination and azimuth, are obtained from considerations involving: (a) a bit/rock interaction law that relates the force and moment acting on the bit to its penetration into the rock; (b) kinematic relationships that describe the local borehole geometry in relation to the bit penetration; and (c) a beam model for the bottom-hole assembly (BHA) that expresses the force and moment at the bit as functions of the external loads applied on the BHA and the geometrical constraints arising from the stabilizers conforming to the borehole geometry. The analytical nature of the propagation equations makes it possible to conduct a systematic stability analysis in terms of a key dimensionless group that controls the directional stability of the drilling system. This group depends on the downhole weight on bit (WOB), on properties of the BHA, on the bit bluntness, and on parameters characterizing the steering response of the bit. The directional stability of a particular system then is assessed by comparing the magnitude of this group to a critical value representing a Hopf bifurcation of stability, which depends only on the BHA configuration and the bit walk. If this group is less than the critical value, the system is referred to as being directionally unstable and borehole spiraling is likely. Stability curves for an idealized BHA with two stabilizers show that the bit walk tends to make drilling systems more prone to spiraling. The influence of the design of push-the-bit rotary steerable systems (RSS) on the onset of spiraling is also discussed, as well as the ability of the stabilizers to tilt freely or not. For directionally unstable systems, the resulting limit cycle, corresponding to the spiral, has been captured by the introduction of the relevant nonlinearity: a saturation of the bit tilt. This nonlinearity enables to characterize further the influence of the directional stability on the borehole geometry and illustrates how the amplitude of the spiral also depends on the borehole inclination and the RSS force. Applications to field cases are discussed and model predictions are tested against actual spiraled-hole data. Simulations conducted by integrating the equations of borehole propagation also are presented. For unstable systems, the model predicts spiraled boreholes with a pitch comparable to what generally is observed in the field. The general good agreement between the model predictions and the field data suggests potential direct implementations. It could lead to model-based control algorithms limiting micro-tortuosity or to better bit and BHA designs, if complemented with field campaigns aiming at refining the model parameters.

Cells defence against stresses that can cause DNA damage (single-strand breaks, double-strand breaks) is crucial in safeguarding the integrity of the genome and the survival of the organism as a whole. One of the genes that plays a pivotal role in maintaining the stability of the genome in humans is p53, which encodes its product p53 protein. The regulation of p53 activation is extremely complex, and molecular cell biology has gathered parts and pieces of the whole pathway. Mental intuition of this complex regulation is challenging; therefore, it requires a different method to quantitatively model and analyse to enhance the current understanding. This thesis has attempted to create two quantitative models of the mechanisms that regulate p53 basal levels and its appropriate activation as a stress response in deciding cell fate by either cell cycle arrest (to stop proliferation of DNA-damaged cells) or apoptosis (programmed cell death) to eliminate damaged cells. In the first part of the research, a modified and improved model from Sun et al. (2011) deterministic model is proposed to explain the p53 basal dynamics and its response to stress due to DNA double-strand breaks. This model in the form of delay differential equations incorporates the most recently found molecular interactions and hypothesis: the core regulators consist of ATM, Mdm2, MdmX, Wip1 and p53. ATM as a stress transducer, amplifies the stress signal and activates p53 and inhibits its regulators Mdm2 and MdmX. The network structure consists of two positive feedback loops (p53 auto-regulation and ATM auto-activation), three negative feedback loops (Mdm2, MdmX and Wip1) and the interplay of p53, Mdm2 and MdmX that have successfully captured the basal dynamics (spontaneous pulses under non-stressed conditions) and stress response (repeated pulses or oscillations under stressed conditions). The model simulation results show that p53 spontaneous pulses are due to intrinsic DNA damage involving low number of DNA double-strand breaks; and p53 auto-regulation is an important positive feedback contributing to a threshold activation of p53 in generating pulses whether spontaneous or repeated. It also shows that p53 dynamics are excitable, in that once initiated, it completes the pulse even if stress signal is inhibited. Bifurcation analysis revealed a spectrum of p53 behaviour under stressed and non-stressed conditions and characterised p53 dynamics as Type II excitability (oscillations arises from non-zero frequency). Most importantly, we reveal some novel findings on the mechanism of threshold activation of p53 pulsatile and oscillatory dynamics that are crucial for its physiological function as a transcription factor and guardian of the genome. The second model is an extension of the first model by incorporating the apoptosis initiation module structure from Zhang et al. (2009a) with modified parameter values for modelling the core regulatory mechanism of p53 protein that activates apoptotic switch in response to high DNA double-strand breaks. The apoptosis initiation module includes Puma, Bcl2 and Bax. p53 activates the transcription of Puma (BH3-only protein that is pro-apoptotic) as a trigger of apoptosis that inhibits Bcl2 protein (pro-survival) and directly activates Bax. Activation of Bax was assumed to be an indicator of apoptosis initiation. The constructed model demonstrated how molecular interactions and stress signal amplification from ATM auto-activation in the p53 network control cell life and death decisions. Particularly, the model simulation results are qualitatively consistent with the experimental findings of an all-or-none activation of apoptosis and predicted overexpression of Bcl2 as a factor in causing the malfunction of the apoptotic switch. This model presents a simplified yet plausible model for molecular mechanism that regulates p53 activation of the apoptotic switch. The model gives insight into the design principles underlying p53 regulation of apoptosis. In summary, the two models presented in this thesis have proposed plausible design principles of p53 basal dynamics and DNA damage response, and activation of apoptotic switch. These models provide novel theoretical insights into p53 regulation.

Bone morphogenic proteins (BMPs) and their inhibitors are increasingly being appreciated for their role in vascular development. Matrix Gla protein (MGP) and crosssveinless 2 (CV2) antagonize BMP-4 and BMP-9, respectively, and the expression of both are stimulated by time-delayed BMP-9/ALK1 signaling, while ALK1 is stimulated by delayed BMP-4 signaling. This double negative feedback system is an important regulator of endothelial cell differentiation, angiogenesis, and maturation. Our collaborators recently discovered oscillatory gene expression of these proteins in human pulmonary artery endothelial cells (HPAECs) undergoing maturation in cell culture. We set out to explore the biological mechanisms guiding this periodic gene expression and to propose a functional role for these oscillations. We used delay differential equations to mathematically model this network and applied dimensionless analysis to study the impact of relative protein degradation rates on oscillations. Monte Carlo parameter sampling and Hopf bifurcation analysis, while adjusting 1 and 2 parameters, helped identify conditions yielding periodic dynamics. We discovered that the presence of the activator, BMP-4, is critical for oscillations and that its inhibitor, MGP, opposes them, making their mutual regulation important for oscillatory gene expression. We also identified other mechanisms for oscillations including the necessary existence of three time-delays, the presence of BMP-9 as a feedback regulator, and the function of CV2 as an indirect inhibitor of MGP. Findings from the cell culture experiments and the mathematical model led us to postulate that oscillations in this network are important for the maintenance of progenitor HPAECs and the lack of oscillations is associated with HPAEC maturation.

A Hamiltonian model for the radial and translational dynamics of clusters of coupled bubbles in an incompressible liquid developed by Ilinskii, Hamilton, and Zabolotskaya [J. Acoust. Soc. Am. 121, 786-795 (2007)] is extended to included the effects of compressibility in the host liquid. The bubbles are assumed to remain spherical and translation is allowed. The two principal effects of liquid compressibility are time delay in bubble interaction due to the finite sound speed and radiation damping due to energy lost to acoustic radiation. The incorporation of time delays produces a system of delay differential equations of motion instead of the system of ordinary differential equations in models of bubble interaction in an incompressible medium. The form of the Hamiltonian equations of motion is significantly different from the commonly used models based on Rayleigh-Plesset equations for coupled bubble dynamics, and it provides certain advantages in numerical integration of the time-delayed equations of motion. Corrections for radiation damping in clusters of interacting bubbles are developed in the form of a time-delayed expression for bubble self-action following the method of Ilinskii and Zabolotskaya [J. Acoust. Soc. Am. 92, 2837-2841 (1992)]. A set of approximate series expansions of this delayed expression is calculated to first order in the ratio of bubble radius to the characteristic wavelength of acoustic radiation from the bubble, and to varying orders in the ratio of bubble radius to characteristic bubble separation distance. Stability of the delay differential equations of motion is analyzed with four successive levels of approximation for the effects of radiation damping and time delay. The stability is analyzed with and without the effects of viscous and thermal damping. The effect of time delay and radiation damping on the pressure radiated by small systems of bubbles is considered. An approximate method to account for the delays in bubble interaction in a weakly compressible liquid is presented. This method converts the system of delay differential equations into an approximate system of ordinary differential equations, which may simplify numerical integration. Several sets of model equations incorporating propagation time delay in bubble interactions are solved numerically with existing algorithms specialized for delay differential equations. Numerical simulations of the dynamics of single bubbles, pairs of bubbles, and clusters of bubbles are used to compare the different levels of approximation for compressibility effects for low- and high-amplitude radial motion in systems of bubbles under free response and pulsed excitation by an external pressure source.

In this thesis we construct a new model of the immune response to the growth of breast cancer cells and investigate the impact of certain drug therapies on the cancer. We use delay differential equations to model the interaction of breast cancer cells with the immune system. The new model is constructed by combining two previous models. The first model accounts for different cell cycles and includes terms to evaluate drug treatments, but ignores quiescent tumor cells. The second model includes quiescent cells, but ignores the immune response and drug treatments. The new model is obtained by combining and modifying these two models to account for quiescent cells, immune cells and includes drug intervention terms. This new model is used to evaluate the effects of pulsed applications of the drug Paclitaxel for models with and without quiescent cells. We use sensitivity equation methods to analyze the sensitivity of the model with respect to the initial number of immune cytotoxic T-cells. Numerical experiments are conducted to compare the model predictions to observed behavior.

In this dissertation, delay differential equation models from mathematical biology are studied, focusing on population ecology. In order to even begin a study of such models, one must be able to determine the linear stability of their steady states, a task made more difficult by their infinite dimensional nature. In Chapter 2, I have developed a method of reducing such questions to the problem of determining the existence or otherwise of positive real roots of a real polynomial. The method of Sturm sequences is then used to make this determination. In particular, I developed general necessary and sufficient conditions for the existence of delay-induced instability in systems of two or three first order delay differential equations. These conditions depend only on the parameters of the system, and can be easily checked, avoiding the necessity of simulations in these cases. With this tool in hand, I begin studying delay differential equations for single species, extending previously obtained results about the existence of periodic solutions, and developing a proof for a previously unproven case. Due to the infinite dimensional nature of these equations, it is quite difficult to prove the existence of periodic solutions. Nonetheless, knowledge of their existence is essential if one is to make decisions about the suitability of such models to biological situations. Furthermore, I explore the effect of delay-dependent parameters in these models, a feature whose use is becoming more common in the mathematical biology literature. Finally, I look at a delayed predator-prey model with delay dependent parameters. Although I was unable to obtain a complete proof for the existence of periodic solutions, significant progress has been made in understanding the nature of this system, and it is hoped that future work will continue to clarify this picture. This model seems to display chaotic behavior for certain parameter regimes, and thus the existence of periodic solutions may be precluded in the most general case.

Genistein is an endocrine-active compound found naturally in soy products. It has been linked to various health effects, both beneficial and adverse. The liver is a major site of genistein transformation. Experimental data suggest genistein is metabolized in the liver into its glucuronide form, and taken back into the gut lumen via biliary excretion. The data show a nontrivial, dose dependent delay in biliary excretion of genistein. Traditional physiologically-based pharmacokinetic (PBPK) modeling methods fail to accurately describe the observed data. We have developed several models that incorporate techniques not typically found in PBPK models to simulate the observed dynamics. The first of these models developed is based on delay differential equations (DDE), where the observed lag in biliary excretion is mathematically described by a time delay. Existence and uniqueness of a solution to the system of equations was obtained and the unknown parameters were obtained via an inverse problem formulation. The nonlinear system of state-dependent delay equations was approximated by a three-stage, implicit, Runge-Kutta method. Using a statistical hypothesis test, we showed that the delay model with the optimal set of parameters obtained is a statistically significant improvement over the PBPK model in simulating the experimental data. Our second modeling approach was taken by considering a dispersion modeling technique. We seeked to develop a fully functioning model for the liver that simulates the distribution, metabolism and excretion of chemicals, and is able to accommodate spatial variations in biologically- and physiologically-based parameters. Our modeling strategy considers the liver as a series of cylindrical tubes, one representing the blood vessel space and one representing the bile duct space, with hepatocytes in between. Dispersion coefficients were adjusted to create a biologically relevant distribution of the concentration of genistein and its metabolites. We established existence and uniqueness of a solution to the system of equations and continuous dependence of the solution on the initial data. A numerical code based on the finite element method was developed to solve the mixed system of nonlinear ordinary and partial differential equations. Parameter estimations were obtained via an inverse problem formulation.