Mathematical modeling with applications in biological systems, physiology, and neuroscience

Doctor of PhilosophyDepartment of MathematicsBacim AlaliDynamical systems modeling is used to describe different biological and physical systems as well as to predict the interactions between multiple components of a system over time. A dynamical system describes the evolution of a given system over time using a set of mathematical laws, typically described by differential equations. There are two main methods to model the dynamical behaviors of a system: continuous time modeling and discrete-time modeling. When the time between two measurements is negligible, the continuous time modeling governs the evolution of the system, however, when there is a gap between any two consecutive measurements, discrete-time system modeling comes into play. Differential equations are used to model continuous systems and iterated maps represent the generations in discrete-time systems. In this dissertation, we study some dynamical systems and present their applications to different problems in biological systems, physiology, and neuroscience. In chapter one, we study the local dynamics of some interesting systems and show the local stable behavior of the system around its critical points. Moreover, we investigate the local dynamical behavior of different systems including the Hénon-Heiles system, the Duffing oscillator, and the Van der Pol equation. Furthermore, we discuss about the chaotic behavior of Hamiltonian systems using two different and new examples. In chapter two, we consider some models in computational neuroscience. Due to the complexity of nerve systems, linear modeling methods are not sufficient to understand the various phenomena in neuroscience. We use nonlinear methods and models, which aim at capturing certain properties of the neurons and their complex dynamics. Specifically, we explore the interesting phenomenon of firing spikes and complex dynamics of the Morris-Lecar model. We consider a set of parameters such that the model exhibits a wide range of phenomenon. We investigate the influences of injected current and temperature on the spiking dynamics of Morris-Lecar model. In addition, we study bifurcations, and computational properties of this neuron model. Moreover, we provide a bound for the membrane potential and a certain voltage value or threshold for firing the spikes. Studying the two co-dimension bifurcations demonstrates more complicated behaviors for this single neuron model. Furthermore, we describe the phenomenon of neural bursting and investigate the dynamics of Morris-Lecar model as a square-wave burster, elliptic burster and parabolic burster. Pharmacokinetic models are mathematical models, which provide insights into the interaction of chemicals with certain biological processes. In chapter three, we consider the process of drug and nanoparticle (NPs) distribution throughout the body. We use a tricompartmental model to study the perfusion of NPs in tissues and a six-compartmental model to study drug distribution in different body organs. We perform global sensitivity analysis by LHS Monte Carlo method using Partial Rank Correlation Coefficient (PRCC). We identify the key parameters that contribute most significantly to the absorption and distribution of drugs and NPs in different organs in the body. In chapter four, we study two infectious disease models and use nonlinear optimization and optimal control theory to help in identifying strategies for transmission control and forecasting the spread of infectious diseases. We analyze the effect of vaccination on the disease transmission in these models. Moreover, we perform global sensitivity analysis to investigate the key parameters in these models. In chapter five, we investigate the complex dynamics of two-species Ricker-type discrete-time competitive model. We perform local stability analysis for the fixed points of the system and discuss about its persistence for boundary fixed points. This system inherits properties of the dynamics of a one-dimensional Ricker model such as the cascade of period-doubling bifurcation, periodic windows, and chaos. We explore the existence of chaos for the equilibrium points for a specific case of this system using Marotto theorem and show the existence of snap-back repeller. In chapter six, we study the problem of chaos synchronization in certain discrete-time dynamical systems. We introduce a drive-response discrete-time dynamical system, which is coupled using convex link function. We investigate a synchronization threshold, after which, the drive-response system uncouples and loses its synchronized behaviors. We apply this method to the synchronized cycles of the Ricker model and show that this model displays a rich cascade of complex dynamics from a stable fixed point and cascade of period-doubling bifurcation to chaos. We numerically verify the effectiveness of the proposed scheme and demonstrate how this type of coupling affects the synchronization of the system. In chapter seven, we study the synchronized cycles of a generalized Nicholson-Bailey model. This model demonstrates a rich cascade of complex dynamics from a stable fixed point to periodic orbits, quasi periodic orbits and chaos. We introduce a coupling of these two chaotic systems with different initial conditions and show how they synchronize over a short time. We investigate the qualitative behavior of Generalized Nicholson-Bailey model and its synchronized model using time series analysis and its long-time dynamics by using its bifurcation diagram

CANCER ASSOCIATED FIBROBLAST DERIVED ANGIOGENIC FACTOR MFAP5 IN OVARIAN CANCER PROGRESSION

Advanced stage ovarian cancer is the most lethal gynecologic malignancy. No major improvement on patient survival has been achieved in the past decade. Therefore, identification of predictive or prognostic markers and further understanding of the molecular mechanisms in ovarian cancer progression are of paramount importance. While cancer cells have always been the targets for the identification of prognostic and predictive markers, the potential for developing new diagnosis and treatments based on the tumor supporting stromal microenvironment is relatively unexplored. Using transcriptome profiling analysis on microdissected stromal and epithelial components of normal and malignant ovarian tissues, we identified a gene signature in the fibroblastic stromal component of the tumor tissue. Among the differentially expressed genes identified in microdissected cancer associated fibroblasts, microfibrillar associated protein 5 (MFAP5), which showed 10 folds increase in expression showed significant association with poorer overall survival in patients with high grade serous ovarian cancer. In addition, we identified a positive correlation between stromal MFAP5 expression and intratumor microvessel density, suggested a pro-angiogenic role of MFAP5. This study aims at delineating the functional roles and mechanisms by which stromal MFAP5 modulates ovarian tumor progression and angiogenesis and investigating the potential of targeting stromal MFAP5 as therapy for ovarian cancer. Functional studies demonstrated that MFAP5 stimulated motility and invasion potential of ovarian cancer cells and microvascular endothelial cells through the binding of αvβ3 integrin receptor. In vivo, targeting stromal MFAP5 using siRNA incorporated chitosan nanoparticles significantly reduced tumor growth, metastasis and intratumoral microvessel density. Further, pathway analyses and western blot analyses demonstrated that CAF-derived MFAP5 modulates ovarian cancer cell motility and invasion potential through the calcium dependent FAK/CREB/TNNC1 signaling pathway and MFAP5 enhanced endothelial cell motility potential and permeability via focal adhesion and stress fiber formation by activating the calcium dependent FAK/ERK/LPP signaling pathway. Finally, monoclonal anti-MFAP5 antibodies were developed. These antibodies demonstrated inhibitory effects on tumor growth and improved paclitaxel bioavailability via the reduction of intratumoral microvessel leakiness. To conclude, our results demonstrated the important roles of MFAP5 in promoting ovarian tumor progression and the potential of targeting stromal MFAP5 as a novel therapeutic approach for ovarian cancer

First-order statistical speckle models improve robustness and reproducibility of contrast-enhanced ultrasound perfusion estimates

Contrast-enhanced ultrasound (CEUS) permits the quantification and monitoring of adaptive tumor responses in the face of anti-angiogenic treatment, with the goal of informing targeted therapy. However, conventional CEUS image analysis relies on mean signal intensity as an estimate of tracer concentration in indicator-dilution modeling. This discounts additional information that may be available from the first-order speckle statistics in a CEUS image. Heterogeneous vascular networks, typical of tumor-induced angiogenesis, lead to heterogeneous contrast enhancement of the imaged tumor cross-section. To address this, a linear (B-mode) processing approach was developed to quantify the change in the first-order speckle statistics of B-mode cine loops due to the incursion of microbubbles. The technique, named the EDoF (effective degrees of freedom) method, was developed on tumor bearing mice (MDA-MB-231LN mammary fat pad inoculation) and evaluated using nonlinear (two-pulse amplitude modulated) contrast microbubble-specific images. To improve the potential clinical applicability of the technique, a second-generation compound probability density function for the statistics of two-pulse amplitude modulated contrast-enhanced ultrasound images was developed. The compound technique was tested in an antiangiogenic drug trial (bevacizumab) on tumor bearing mice (MDA-MB-231LN), and evaluated with gold-standard histology and contrast-enhanced X-ray computed tomography. The compound statistical model could more accurately discriminate anti-VEGF treated tumors from untreated tumors than conventional CEUS image. The technique was then applied to a rapid patient-derived xenograft (PDX) model of renal cell carcinoma (RCC) in the chorioallantoic membrane (CAM) of chicken embryos. The ultimate goal of the PDX model is to screen RCC patients for de novo sunitinib resistance. The analysis of the first-order speckle statistics of contrast-enhanced ultrasound cine loops provides more robust and reproducible estimates of tumor blood perfusion than conventional image analysis. Theoretically this form of analysis could quantify perfusion heterogeneity and provide estimates of vascular fractal dimension, but further work is required to determine what physiological features influence these measures. Treatment sensitivity matrices, which combine vascular measures from CEUS and power Doppler, may be suitable for screening of de novo sunitinib resistance in patients diagnosed with renal cell carcinoma. Further studies are required to assess whether this protocol can be predictive of patient outcome

Modification of left ventricular geometry and function during healing after acute myocardial infarction

Increased left ventricular (LV) size and deformation of LV geometry are associated with LV dysfunction. Regional shape distortion (RSD), detected on two-dimensional echocardiography (2D-Echo) after acute myocardial infarction (MI), is associated with poor outcome. Two hypotheses were tested: i) early RSD of the asynergic infarct zone after MI is followed by progressive global LV dilatation, remodelling towards a spheroidal shape, and more LV dysfunction; and ii) the progressive remodelling of LV geometry spans the phases of early infarction and healing and may be modified by early and prolonged therapies applied over the phases of infarction and healing. A bench to bedside approach was used, with concurrent studies in a dog model of healing over 6 weeks after MI and patients with first MI's. Computer- assisted analysis of the 2D-Echo images with 3D reconstruction was used to quantify LV asynergy (akinesis + dyskinesis), LV volumes, LV ejection fraction, RSD bulge and global LV shape. The animal studies showed that collagen deposition during healing after MI increases progressively, reaching a plateau around 2 weeks, and deposition of collagen in already dilated infarct zones is followed by late thinning and further RSD associated with LV aneurysms. Importantly, serial 2D-Echo tracked the in- vivo changes in LV geometry and function and showed greater RSD and LV dysfunction with anterior than inferior MI, and with transmural MI than nontransmural MI. Other studies showed: i) lower LV resistance to distension and rupture in infarcted hearts; ii) marked extracellular matrix (ECM) disruption and RSD in transmural MI; ill) delayed effects on LV remodelling after infarct-limiting therapies given during acute MI; iv) loss of beneficial effects of the vasodilator nitroglycerin (NTG) with hypotension induced by high doses during acute MI; v) decreased wall stress by prolonged LV unloading after MI, with nitrates (eccentric dosing) and angiotensin-converting enzyme (ACE) inhibitors, limited early RSD and progressive LV remodelling and dysfunction; this effect was greater with therapy over 6-weeks than just over the first 2 weeks; vi) late reperfusion limited early RSD and adverse LV remodelling, and preserved ECM in the epicardial rim; vii) the resistance of the healed left ventricle to distension and rupture was further reduced by prolonged anti-inflammatory therapy (ibuprofen); viii) prolonged ACE inhibitor therapy decreases infarct collagen, which may be harmful under certain conditions. The clinical studies with serial 2D-Echo showed that systematic tomographic imaging could provide quantitative data on regional and global LV geometry and function including the degree of RSD (depth, area, and volume). An early 2D-Echo not only provided diagnostic data on LV thrombi and complications of MI, but the extent of LV asynergy on the initial 2D-Echo predicted outcome at 3 months and 1 year. Importantly, the degree of RSD on the initial 2D-Echo predicted patients at high risk of adverse remodelling with infarct expansion, greater LV dysfunction, progressive LV dilatation, and poor outcome at 1 year. Survivors of MI with > 18% LV asynergy and significant RSD on a baseline 2D Echo were at increased risk of topographic deterioration on exercise programs. Anti-inflammatory therapy after MI resulted in more RSD and adverse remodelling. Short-term LV unloading with low-dose intravenous NTG therapy during the acute MI, as well as prolonged nitrate (eccentric dosing) and captopril therapy during healing over 6 weeks after MI, improved 2D-Echo indexes of LV geometry and function, decreased complications and improved outcome. Acute thrombolytic therapy also limited LV remodelling after MI. In all these studies, the degree of RSD and severity of LV dysfunction were greater with anterior than inferior MI, and with Q-wave than non-Q wave MI. In Conclusion, the overall results indicate that early RSD in the infarct zone leads to progressive global LV dilatation, LV dysfunction and poor outcome and the changes in LV geometry and function can be quantified by serial quantitative 2D-Echo imaging. Marked RSD is associated with early ECM disruption and aneurysm formation after transmural MI. During healing, infarct zones may be thinned and dilated before the collagen plateau, and collagen deposition into these zones result in further RSD and chronic aneurysms. Prolonged anti-remodelling therapy during healing, with agents that decrease wall stress without damaging the ECM, or decreasing infarct collagen, or causing infarct thinning, or impairing healing, might be more effective for reducing RSD, LV aneurysm, global dilatation and poor outcome. The 2D-Echo measurement of RSD early after MI might be potentially important for stratifying patients according to their topographic status and for the objective assessment of the effects of anti-remodelling strategies during healing after MI

Mechanistic insights into neuronal oscillatory activity in the dopamine-intact and dopamine-depleted primary motor cortex

In Parkinson’s disease (PD) the loss of the neurotransmitter dopamine (DA) results in abnormal oscillations of the cortico-basal ganglia network, the emergence of which correlate with symptoms. Increased oscillatory power in the primary motor cortex (M1) is reduced by dopamine replacement therapy and by targeted stimulation, suggesting that M1 plays an important role in the pathology of PD. In this study we have investigated, using pharmacology, the mechanisms by which oscillatory activity in rat M1 is generated and determined the power changes associated with DA depletion and DA receptor modulation. Extracellular local field potential recordings were made in brain slices of M1 which were prepared using a modified protocol to improve viability. Co-application of carbachol (5 μM) and kainic acid (100 nM) elicited simultaneous theta (4-8 Hz) and gamma (30-40 Hz) oscillations in layer V of M1. These oscillations displayed phase-amplitude coupling; the first report of such findings in vitro. These oscillations were found to be pharmacologically distinct with theta oscillations generated by intrinsic non-synaptic mechanisms while gamma oscillations required contributing excitatory and inhibitory networks. Following successful unilateral lesions using 6-hydroxydopamine (6-OHDA), as determined by the adjusting step test, DA-depleted (ipsilateral) and DA-intact (contralateral) slices were obtained. Although no difference in the oscillatory profile of M1 ipsilateral, contralateral or age-matched control (AMC) slices was found, bath application of DA reduced gamma power only in the ipsilateral slices and amphetamine only decreased gamma power in contralateral slices. Furthermore, D2-like receptor activation consistently increased both theta and gamma power in contralateral and AMC slices, while only theta power was increased in ipsilateral slices. Overall, these data suggest that DA, through action at multiple sites, differentially modulates the power of both theta and gamma oscillations in M1. Using the 6-OHDA model, the oscillatory activity of M1 in vivo was investigated. Successful lesions were determined by using the rotometer, the locomotor activity and the adjusting stepping tests at 2-4 weeks post-surgery. Further testing at 22 weeks post-surgery indicated the long-term stability of the lesions. Using depth electrode and EEG recordings, oscillatory activity in the 2-10 Hz range was found in the ipsilateral and contralateral hemispheres of both lesioned and sham animals. However, only in the ipsilateral hemisphere of DA-depleted animals did we detect a 30-40 Hz oscillatory peak, which was localised to layer V of M1. In EEG recordings this led to a significant increase in the interhemispheric ratio. Using depth electrode recordings, the ipsilateral 30-40 Hz oscillation (but not 2-10 Hz oscillation) was reduced by the administration of L-DOPA (6 mg/kg) with a reduction in interhemispheric ratio. However, administration of zolpidem (0.3 mg/kg), which previously reduced abnormal beta oscillatory activity in vivo and in vitro resulting in the rebalancing of interhemispheric beta power (Hall et al., 2014; Prokic et al., 2015), was without effect. Overall, these studies demonstrate that M1 alone can generate multiple, pharmacologically distinct, but interacting oscillations, which contribute to pathological activity in the DA-depleted state

2011 IMSAloquium, Student Investigation Showcase

Inquiry Without Boundaries reflects our students’ infinite possibilities to explore their unique passions, develop new interests, and collaborate with experts around the globe.https://digitalcommons.imsa.edu/archives_sir/1003/thumbnail.jp

Modernizing Markov Chains Monte Carlo for Scientific and Bayesian Modeling

The advent of probabilistic programming languages has galvanized scientists to write increasingly diverse models to analyze data. Probabilistic models use a joint distribution over observed and latent variables to describe at once elaborate scientific theories, non-trivial measurement procedures, information from previous studies, and more. To effectively deploy these models in a data analysis, we need inference procedures which are reliable, flexible, and fast. In a Bayesian analysis, inference boils down to estimating the expectation values and quantiles of the unnormalized posterior distribution. This estimation problem also arises in the study of non-Bayesian probabilistic models, a prominent example being the Ising model of Statistical Physics. Markov chains Monte Carlo (MCMC) algorithms provide a general-purpose sampling method which can be used to construct sample estimators of moments and quantiles. Despite MCMC’s compelling theory and empirical success, many models continue to frustrate MCMC, as well as other inference strategies, effectively limiting our ability to use these models in a data analysis. These challenges motivate new developments in MCMC. The term “modernize” in the title refers to the deployment of methods which have revolutionized Computational Statistics and Machine Learning in the past decade, including: (i) hardware accelerators to support massive parallelization, (ii) approximate inference based on tractable densities, (iii) high-performance automatic differentiation and (iv) continuous relaxations of discrete systems. The growing availability of hardware accelerators such as GPUs has in the past years motivated a general MCMC strategy, whereby we run many chains in parallel with a short sampling phase, rather than a few chains with a long sampling phase. Unfortunately existing convergence diagnostics are not designed for the “many short chains” regime. This is notably the case of the popular R statistics which claims convergence only if the effective sample size per chain is large. We present the nested R, denoted nR, a generalization of R which does not conflate short chains and poor mixing, and offers a useful diagnostic provided we run enough chains and meet certain initialization conditions. Combined with nR the short chain regime presents us with the opportunity to identify optimal lengths for the warmup and sampling phases, as well as the optimal number of chains; tuning parameters of MCMC which are otherwise chosen using heuristics or trial-and-error. We next focus on semi-specialized algorithms for latent Gaussian models, arguably the most widely used of class of hierarchical models. It is well understood that MCMC often struggles with the geometry of the posterior distribution generated by these models. Using a Laplace approximation, we marginalize out the latent Gaussian variables and then integrate the remaining parameters with Hamiltonian Monte Carlo (HMC), a gradient-based MCMC. This approach combines MCMC and a distributional approximation, and offers a useful alternative to pure MCMC or pure approximation methods such as Variational Inference. We compare the three paradigms across a range of general linear models, which admit a sophisticated prior, i.e. a Gaussian process and a Horseshoe prior. To implement our scheme efficiently, we derive a novel automatic differentiation method called the adjoint-differentiated Laplace approximation. This differentiation algorithm propagates the minimal information needed to construct the gradient of the approximate marginal likelihood, and yields a scalable differentiation method that is orders of magnitude faster than state of the art differentiation for high-dimensional hyperparameters. We next discuss the application of our algorithm to models with an unconventional likelihood, going beyond the classical setting of general linear models. This necessitates a non-trivial generalization of the adjoint-differentiated Laplace approximation, which we implement using higher-order adjoint methods. The generalization works out to be both more general and more efficient. We apply the resulting method to an unconventional latent Gaussian model, identifying promising features and highlighting persistent challenges. The final chapter of this dissertation focuses on a specific but rich problem: the Ising model of Statistical Physics, and its generalization as the Potts and Spin Glass models. These models are challenging because they are discrete, precluding the immediate use of gradient-based algorithms, and exhibit multiple modes, notably at cold temperatures. We propose a new class of MCMC algorithms to draw samples from Potts models by augmenting the target space with a carefully constructed auxiliary Gaussian variable. In contrast to existing methods of a similar flavor, our algorithm can take advantage of the low-rank structure of the coupling matrix and scales linearly with the number of states in a Potts model. The method is applied to a broad range of coupling and temperature regimes and compared to several sampling methods, allowing us to paint a nuanced algorithmic landscape