The sine-Gordon model was modified to depict the dynamics of the double helix of a DNA system. The model was based on an the inhomogeneities that exist in the base sequence of the DNA structure. Various works that has been done in similar areas were discussed as well as the method that was used . Computer simulations were performed on the model to see the distortions that occurs in the DNA system
The coordinated migration of multiple cells in a tissue is critical to a variety of biological processes, such as wound healing, cancer invasion, and morphogenesis. During vertebrate gastrulation, a migratory population of mesoderm and endoderm, collectively referred to as mesendoderm, migrates across the fibronectin-rich blastocoel roof (BCR) of the embryo. Leading-edge mesendoderm cells exhibit polarized protrusive behaviors, integrin-mediated tractions and directional migration along the BCR. Leading row traction forces are balanced by c-cadherin dependent cell-cell adhesions, which are required to pull follower row cells forward. Studies of collective cell migration in gastrulation have focused in large part on tissue explants removed from Xenopus laevis embryos. However, some cell behaviors (e.g., migration speed) change when the tissues are removed from the embryo and placed in vitro, confounding efforts to understand mechanisms of cell migration and tissue formation. To help address these limitations we are developing in silico approaches. We have constructed an agent based model (ABM) in the Cellular-Potts framework to investigate collective cell migration in the Xenopus embryo. Our model consists of 9 rules governing leader and follower cell dynamics and represents a biological dorsal marginal zone (DMZ) explant of 64 cells with both leader and follower cell agents over a 2 hour timecourse. Our model was calibrated to reproduce published experiments demonstrating mechanical properties of anisotropic tension in the DMZ explant. Model predictions suggests that cell intercalation and in vivo geometry contribute to increased collective cell migration speed during mesendoderm mantle closure along the BCR during gastrulation.
Vision play vital roles in the lives of every animal. Rods and cones are two primary photoreceptor cells in the eye responsible for converting light energy (photon) into electrical signal perceived by the brain to enable vision. To prevent accumulation of toxics caused by photo-oxidative compounds, the rod and cone cells undergo daily renewal through addition of new disks at the base of their outer segment and removal of older ones from the tip. The balance between these two processes help the cell to maintain constant or an equilibrium length necessary for optimal performance of these cells. Imbalance may lead to retinal disease such as retinitis pigmentosa, a form of inherited blinding disease caused by degeneration of rod cells followed by progressive lost of cone cells. Also, when the retina is detached from the retinal pigmented epithelium (RPE), the rod and cone cells degenerate and if the retina is reattached on time, the cells are able to regenerate to restore vision. When the rod outer segment suffer from degeneration due to retinal detachment, at what point will regeneration be impossible? How does retinal detachment disrupt renewal process (addition of new disks and shedding)? What mechanism controls the renewal process? How can degenerating rod and cone cells be rescue? These are some of the questions we intend to quantitatively address using mathematical model and comparing the result to a date obtained from zebra-fish. We focused on rod cells because survival of cone cells depends on rod cells and also the disks in the rod outer segment are discrete except few newly formed disks at the base which are still connected to the cell membrane and to one another. This feature make rod cells easily trackable and obtaining experimental data quiet easier compare to cone cells.
: Jurkat cells are an immortalized line of human T lymphocytes that are commonly used to study leukemia, HIV, and Calcium (Ca2+) signaling. Stimulation of Jurkat cells leads to the depletion of the cells internal storage of Ca2+, the Endoplasmic Reticulum (ER), which in turn causes Ca2+ to enter the cell via a mechanism of Store Operated Calcium Entry (SOCE). The combination of these mechanisms causes oscillations in the Ca2+ concentration of the cell. New data obtained last year, however, suggests that Jurkat cells, in addition to presenting SOCE-based oscillations, are also capable of producing ER-based oscillations when the external entry mechanism is genetically eliminated. With this in mind, I construct a model of ordinary differential equations that incorporates both oscillatory mechanisms and can produce both ER and SOCE based oscillations. ER-based oscillations have been studied before in other cell types, but the interplay between the two oscillatory mechanisms is not well understood. In my poster I present the construction of the model, as well as the defining characteristics of each oscillation; finally, I present a hypothesis for how the interaction between the oscillatory mechanisms results in a dominating signal being exhibited under different conditions.
Parasites exhibit remarkable diversity in their life history traits to adapt to the unique ecological challenges posed by their hosts. Within the genus Plasmodium, the life cycle of the malaria-causing species involves multiple rounds of replication, with a fraction of infected red blood cells being committed to producing specialized stages for onward transmission to vectors. The rate of proliferation is limited by the burst size or the average number of daughter cells to emerge from each infected red blood cell. As proliferation is crucial for establishing and maintaining the infection, parasites would be expected to evolve to the maximal burst size that does not prematurely end the infection by killing its host. In reality, observed burst sizes vary significantly across species and even among strains, suggesting that maximizing the burst size is not always the best strategy. More specifically, restricting within-host proliferation may be beneficial for the parasites though the exact mechanism is unclear.
Using a within-host model parameterized for the rodent malaria, Plasmodium chabaudi, we investigate how host mortality and resource limitation affect the optimal burst size. We focus on the acute phase which encompasses the first and typically largest wave of parasite abundance with most of the parasite’s transmission success gained disproportionately in this phase. By calculating the cumulative transmission potential at the end of the acute phase, we find that the most transmissible strain does not maximize its burst size even if the value does not induce host mortality . Greater proliferation leads to the production of more sexual forms, but there are diminishing returns in transmission success. Moreover, the benefits of faster proliferation come at the cost of significantly shortening the period of high infectivity. Therefore, the optimal burst size emerges from the trade-off between the length of the acute phase and the production of the sexual forms. By identifying resource availability as a key mechanism limiting the burst size, we are better able to understand how parasite traits can influence the varying virulence we see in malaria infections.
Subaerial biofilms (SABs) are well-organized self-sufficient communities that colonize stone surfaces exposed to the atmosphere. Such biofilms are composed of different microbial species embedded in a self-produced matrix of extracellular polymeric substances (EPS), which spreads onto the substratum contributing to the microorganisms protection from external factors. Microbial life within these ecosystems is hard and mainly depends on the availability of liquid water, which plays an essential role in the microbial metabolic activities. Understanding the relationship between ambient air, biofilm and stone is of paramount importance for both stone conservation and biofilm lifecycle. In this talk, a mathematical model describing the water transport through atmosphere, SAB and substratum is presented, taking into account also the effect of the water content of the three layers on the metabolic activities of the microbial communities constituting the biofilm. Numerical simulations are performed to explore how SAB affects the flux of water between atmosphere and substratum. Simulation results are presented and discussed.
Stochastic fluctuations are ubiquitous in natural and man-made systems. Those fluctuations can give rise to dramatic, unexpected, and oftentimes catastrophic dynamical consequences such as a sudden population collapse to extinction. Those events are very rare and never happen on a realistic time scale. We are keenly interested in speeding up such a fluctuation-induced rare event of population extinction. We consider a stochastic Verhulst population growth model in a temporally modulated extrinsic condition. In the absence of temporal extrinsic perturbation, the stochastic population system transits to an extinction state along the optimal path which is a heteroclinic orbit connecting the extinction state to the fluctuation-induced state which is created purely due to stochastic fluctuations. When the temporally modulated extrinsic perturbation is turned on, the population extinction accelerates with its mean passage time to extinction being exponentially inversely proportional to the amplitude of temporal modulation. However, such an acceleration is limited only to the small amplitude temporal modulation beyond which the optimal path is disconnected, making the fluctuation-induced extinction implausible. We enforce the connectedness of the optimal path during large amplitude temporal modulation and thus maximally accelerate the population extinction, by using the (quantum) counter-diabatic control that is able to drive the quantum system in a finite time while keeping the system in the quasi-equilibrium state and suppressing non-adiabatic transitions. We extend its application to a tumor growth model with therapy-induced resistance.
The presence of multiple timescales in complex biological or ecological systems is ubiquitous in nature. Mathematically, these systems are known as 'slow-fast' systems or singularly perturbed systems. Recently, there has been a rising interest among researchers to study the ecological implications of slow-fast systems. The presence of multiple timescales in a biological system inevitably gives rise to the study of long transients. Here, we consider a slow-fast predator-prey model with Bazykin-type reaction kinetics to study the spatio-temporal long transients. The temporal counterpart of the system shows the existence of peculiar periodic solutions: canard and relaxation oscillation. However, a parametric domain is identified where the system shows the existence of two canard cycles, stable and unstable. Even in the spatially extended system, a spatio-temporal canard explosion is observed. This implies that the system dynamics change abruptly from small to large amplitude oscillations within an exponentially small parameter interval. The system dynamics become much more complex near a bifurcation threshold. We argue that the spatial average of the species density over time is not enough to capture the spatial heterogeneity of the distribution of the species. Hence, we introduce two additional metrics to identify the rich spatio-temporal dynamics, which include a variety of long transient regimes.
Since its invention by two trailblazing biologists in 1943, the Luria-Delbruck experiment has been a preferred tool for measuring microbial mutation rates in the laboratory. Practical algorithms for computing a variety of mutant distributions induced by the Luria-Delbruck experiment play a pivotal role in helping biologists obtain accurate estimates of mutation rates. This presentation focuses on an important type of the Luria-Delbruck distribution that simultaneously accommodates differential fitness between mutants and nonmutants and imperfect plating efficiency.
This distribution was earnestly tackled in the 1990s, and important intermediate results were obtained. However, a workable algorithm remained an unachieved expectation at the time. In the 2010s, a clever contour integration approach was taken. An elegant algorithm relying on numerical integration was then devised. Illustrative testing examples showed remarkable performance of the integration-based algorithm. But real-world research problems can be far more challenging than artificial testing examples, and the integration-based method performed dishearteningly on some real-world examples.
I here present an alternative algorithm that effectively exploits some properties of the hypergeometric function. Reliant on the hypergeometric function and simple arithmetic operations, the new algorithm may appear at first sight to be clumsy but computes the mutant distribution more accurately and efficiently. Examples are given to show the usefulness of the new algorithm in actual microbial mutation research.
Humans provide food for wildlife for conservation and recreational purposes, but the resulting aggregation of wildlife around food sources can increase transmission of infectious diseases. Past work investigating the consequences of wildlife feeding for parasite transmission typically assumes that food is provided at a constant rate, but in reality, observations of changing wildlife abundance or infection can influence how much food is provided, potentially generating feedbacks between human behavior and wildlife disease. Focusing on backyard bird-feeding as a widespread and globally popular form of wildlife feeding, I develop a simple mathematical model for coupling the abundance and infection of birds with the intensity of food provisioning, contrasting scenarios where the rate of food provisioning is independent of, or depends on, components of the natural system. Unlike constant food provisioning, which usually results in a small outbreak followed by a smooth approach to equilibrium infection prevalence, coupling food provisioning to bird abundance and infection can result in more complex emergent dynamics, including larger and more frequent epidemic peaks and higher cumulative bird mortality. Accounting for this coupling of human activity to observations of ecological dynamics could inform development of best practice guidelines for wildlife feeding that minimize its unintended negative consequences.
Wolbachia intracellular bacteria successfully reduce the transmissibility of arthropod-borne viruses (arboviruses) when introduced into virus-carrying vectors such as mosquitoes. Despite the progress made by introducing Wolbachia bacteria into the Aedes aegypti wild-type population to control arboviral infections, reports suggest that heat-induced loss-of-Wolbachia-infection as a result of climate change may reverse these gains. Novel, supplemental Wolbachia strains that are more resilient to increased temperatures may circumvent these concerns, and could potentially act synergistically with existing variants. In this work, we model the ecological dynamics among three distinct mosquito (sub)populations: a wild-type population free of any Wolbachia infection; an invading population infected with a particular Wolbachia strain; and a second invading population infected with a distinct Wolbachia strain from that of the first invader. We explore how the range of possible characteristics of each Wolbachia strain impacts mosquito prevalence. Our results show that releasing mosquitoes with two different strains of Wolbachia did not increase their prevalence, compared with a single-strain Wolbachia-infected mosquito introduction and only delayed Wolbachia dominance.
The gut microbiota is critical for maintaining human health, yet chronic intake of certain medications, such as polyethylene glycol (PEG), has the potential to perturb its equilibrium. In this study, we sought to elucidate the competitive dynamics between two bacterial families, Muribaculaceae and Bacteroidaceae, under the influence of altered osmolality due to protracted PEG consumption. Employing the classical competitive exclusion model, we scrutinized variations in the interaction and growth rates of these bacterial taxa as a function of increasing PEG concentrations. Our findings demonstrate that escalating PEG levels provoke significant alterations in the composition of commensal bacteria, with Muribaculaceae being disproportionately affected. The competitive interplay between the bacterial taxa is predominantly governed by their interaction rate; a heightened interaction rate results in intensified competition, corroborating our hypothesis. Muribaculaceae's elevated interaction rate is posited as the primary factor underlying its observed decline in abundance. This study not only provides a deeper understanding of the mechanisms through which PEG consumption disrupts gut microbiota homeostasis but also paves the way for future investigations focusing on the development of targeted interventions to counteract these detrimental effects.
During evolutionary transitions in individuality, new levels of selection are introduced, and thus, the possibility of discordant selection between levels. Cancer or ‘cheating cells’ comes hand in hand with multicellularity. On a cellular level there is selection for cells that replicate faster, even if it is to the detriment of the organism. Given this phenomenon, what has facilitated and maintained the transition to multicellular life? The unicellular bottleneck (unicellular offspring) has been hypothesized as an adaptation to facilitate cooperation among cells by purging lineages of cheating cells. The evolution of propagule size has been explored previously but here we introduce spatial structure and different modes of dispersal by simulating individuals competing on a lattice. We find that size dependent competition favours binary fragmentation strategies but high mutation rates to cancer cells favours fragmentation modes with more unicellular offspring. Specifically, multiple fission is favoured by global dispersal and the unicellular propagule strategy is favoured by local dispersal. Our simulation sheds light on the evolution of multicellular life cycles and the prevalence of unicellular offspring in multicellular species.
We consider a hybrid model of an annual species with the timing of a stage transition governed by density dependent phenology. We show that the model can produce a strong Allee effect as well as overcompensation. The density dependent probability function that describes how population emergence is spread over time plays an important role in determining population dynamics. Our extensive numerical simulations with a density dependent gamma distribution indicate very rich population dynamics, from stable/unstable equilibria, limit cycles, to chaos.
We consider the effect of network structure on the evolution of a population. Models of this kind typically consider a population of fixed size and distribution. Here we consider eco-evolutionary dynamics where population size and distribution can change through birth, death and migration, all of which are separate processes. This allows complex interaction and migration behaviours that are dependent on competition. For migration, we assume that the response of individuals to competition is governed by tolerance to their group members, such that less tolerant individuals are more likely to move away due to competition. We looked at the success of a mutant in the rare mutation limit for the complete, cycle and star networks. Unlike models with fixed population size and distribution, the distribution of the individuals per site is explicitly modelled by considering the dynamics of the population. This in turn determines the mutant appearance distribution for each network. Where a mutant appears impacts its success as it determines the competition it faces. For low and high migration rates the complete and cycle networks have similar mutant appearance distributions resulting in similar success levels for an invading mutant. A higher migration rate in the star network is detrimental for mutant success because migration results in a crowded central site where a mutant is more likely to appear.
Tuberculosis (TB) caused over 1.6 million deaths in 2021. TB is associated with granulomas, organized structures of immune cells that contain the causative bacteria. These structures are three-dimensional with an inner core of macrophages and an outer cuff of T cells. Advanced 3D cell cultures have been applied to emulate these clinical structures in vitro. One in vitro approach showed that 3D spheroid models have improved bacterial control compared to traditional in vitro infection models. We use hybrid modeling to simulate these spheroid models and traditional counterparts, with an agent-based model of immune cell and bacteria rules coupled to a partial differential equation model of chemokine diffusion. We calibrate our model to experimental data while enforcing shared parameters between the spheroid and traditional setups, only changing the initial structure and movement rules to reflect the experimental setups. Lower bacterial load in spheroid simulations as compared to traditional culture simulations is predicted to be due to increased proportions of activated macrophage killing of bacteria, either in tandem with increased proportions of CD8+ T cell activation or not. The spatial distribution of cells was found to be an important factor in macrophage and CD8+ T cell activation in spheroid simulations, with more activation being associated with increased proximity. Next, an in silico experiment was performed, where the initial structure and movement rules were uncoupled to see if either of these independently lead to bacterial control. Neither of the uncoupled mechanisms reduced bacterial load on its own, rather they worked together synergistically. This work further emphasizes the impacts of spatial organization and dimension in biological processes, while highlighting the flexibility of in silico modeling and the perturbations it makes possible.
Equine Infectious Anemia Virus (EIAV) is an important infection in equids, and its similarity to HIV creates hope for a potential vaccine. We analyze a within-host model of EIAV infection with antibody and cytotoxic T lymphocyte (CTL) responses. In this model, the stability of the biologically relevant endemic equilibrium, characterized by the coexistence of long-term antibody and CTL levels, relies upon a balance between CTL and antibody growth rates, which is needed to ensure persistent CTL levels. We determine the model parameter ranges at which CTL and antibody proliferation rates are simultaneously most influential in leading the system towards coexistence and can be used to derive a mathematical relationship between CTL and antibody production rates to explore the bifurcation curve that leads to coexistence. We employ Latin hypercube sampling and least squares to find the parameter ranges that equally divide the endemic and boundary equilibria. We then examine this relationship numerically via a local sensitivity analysis of the parameters. Our analysis is consistent with previous results showing that an intervention (such as a vaccine) intended to control a persistent viral infection with both immune responses should moderate the antibody response to allow for stimulation of the CTL response. Finally, we show that the CTL production rate can entirely determine the long-term outcome, regardless of the effect of other parameters, and we provide the conditions for this result in terms of the identified ranges for all model parameters.
The 2019 novel coronavirus, SARS-CoV-2, is a pathogen of critical significance to international public health. Knowledge of the interplay between molecular-scale virus-receptor interactions, single-cell viral replication, intracellular-scale viral transport, and emergent tissue-scale viral propagation is limited. Moreover, little is known about immune system-virus-tissue interactions and how these can result in low-level (asymptomatic) infections in some cases and acute respiratory distress syndrome (ARDS) in others, particularly with respect to presentation in different age groups or pre-existing inflammatory risk factors. Given the nonlinear interactions within and among each of these processes, multiscale simulation models can shed light on the emergent dynamics that lead to divergent outcomes, identify actionable “choke points” for pharmacologic interventions, screen potential therapies, and identify potential biomarkers that differentiate patient outcomes. Given the complexity of the problem and the acute need for an actionable model to guide therapy discovery and optimization, we introduce an ABM model of SARS-CoV-2 dynamics in lung tissue in PhysiCell. This model finds key points on infection dynamics identified for both early timescale (interferon and receptor uptake) and later lymph node interactions (CD4+/CD8+ T Cell and anti-body recruitment). Interestingly, the model also pointed towards the spatial ability of immune cell sensing and initial infection spread with similar MOIs. More broadly, this effort created a reusable, modular framework for studying viral replication and immune response in tissues, which can also potentially be adapted to related problems in immunology and immunotherapy.
The therapy of cancer is a long-standing and worldwide issue. Since cancers metastasize to other organs, the treatment method of multiple organs simultaneously is required but it is difficult and has not yet been established. In this study, we investigated an optimal therapy method that targets Notch signaling network which shown in multiple cancers in common. As example studies, we targeted embryonal brain tumor (EBT) and non-small cell lung cancer (NSCSC). Both the two cancers undergo oncogenic development through increased HES-1 via Notch signaling, but their signaling pathways of Notch 1 and Notch 2 to enhance HES-1 gene have contrastive roles. In NSCSC, Notch 1/2 activates/inhibits cancerization of cells. In contrast, Notch 1/2 plays an opposite role in EBT, namely, Notch 1/2 inhibits/activates cancerization of cells. To find a possible therapy by which we can treat both cancers at the same time, we developed a conceptual mathematical model based on Notch signaling with the opposite pathways. We explored which network pathway is critical to enhance the cancer cells by sensitivity analysis and found that an intra-cellular pathway is more critical than inter-cellular pathway in enhancing the cancerization of cells and the pathway of Notch transport pathway from cytosol to membrane can be a common network to enhance the cancerization of cells in both cancers. Based on these observations, we also carried out in silico therapy tests for ten patient cases and found that network enhancement therapy is more effective than network cleavage therapy to reduce the number of cancer cells and multiple network therapies are more effective than a treatment of single network therapy. This study suggests that there are optimal signaling network therapies that can treat multiple cancers with contrasting Notch networks and that the simultaneous use of drugs that regulate multiple signaling networks may be most effective in reducing cancer cells.
During the COVID-19 a pandemic, mathematicians mobilized to create models to predict the rise of COVID-19 through communities. In parallel to the spread of the virus, there was an equally insidious spread of misinformation across various social media platforms. In this poster, we will analyze the similarities and differences in transmission of various types of COVID-19 misinformation (e.g, conspiratorial and non-conspiratorial) via semi-viral tweets in the early stages of the pandemic. We build and analyze follower/followee network graphs for retweets of different types of misinformation and determine the characteristics that distinguish the spread conspiratorial versus non-conspiratorial misinformation.
In Bogotá, Colombia, more than 1.8 million cases of COVID-19 and 30,000 deaths had been reported by April 2023. During the critical phase of the pandemic, decision makers required estimates of the impact of different scenarios to design public-health interventions, such as school closures, face-masks, or the distribution of available vaccines. For instance, public schools were closed for in-person instruction in Bogotá during most of 2020. We used an agent-based model of COVID-19 and calibrated it to represent the epidemiological dynamics of COVID-19 in Bogotá, including SARS-CoV-2 variants, and capable of reproducing time-varying public health interventions, such as reduction in mobility, school closures, and vaccination programs. To inform school reopening during the first semester of 2021, we simulated school reopening at different capacities, and found that school reopening could have had a small impact on the number of deaths reported in the city during the third wave at 35% capacity of in-person instruction during the simulation period. Deaths were lowest when only reopening pre-kinder grades, and largest when secondary school was opened. The impact of opening pre-kinder at 100% capacity was below 10% of additional deaths. Finally, we also estimated the impact of vaccination in the city during the third and largest wave of COVID-19 in 2021. Our simulation results suggest that vaccination may have prevented more than 17 thousand deaths in the city.
In this study, we aim to investigate the nature of bifurcations in an extended version of an SVEIR type compartmental model with differential morbidity. Since all existing COVID-19 vaccines are imperfect, we consider vaccine efficacy as a pivotal parameter in the study. The endemic equilibrium of the model was analyzed by explicitly constructing a quadratic equation which was then manipulated appropriately in order to derive R0 using an alternative approach. This also permitted a comprehensive categorization of the number of endemic equilibria based on the threshold condition R0 = 1, which also seemed to imply potential existence of the backward bifurcation phenomenon. However, numerical simulations and application of center manifold theory showed that the bifurcation at R0 = 1 is forward. Therefore, an analysis based on the existence of a quadratic equation at the endemic equilibria is not sufficient in establishing backward bifurcations. We then explored what causes the equation of endemic equilibria to become quadratic and the bifurcation diagram to have a non-linear shape. In this respect, it was shown that the underlying equation is not quadratic (but linear) when the vaccine is perfect which also yields a linear bifurcation diagram.
Keywords: COVID-19 vaccination, Compartmental models, Basic reproduction number, Quadratic equation of endemic equilibria, Bifurcations
During emergent outbreaks of viral infections, public health policy decisions are made on the basis of incomplete information in a changing landscape of scientific knowledge and budgetary and infrastructure constraints. Accounting for the trade-offs necessitated by the resource limitation is essential when formulating an optimal policy response. In this work, we pose optimal control problems to explore the implications of several such trade-off, focusing on testing vs. vaccination and long-term vs. short-term public health objectives. We also explore the how these optimal controls are influenced by the efficacy of the interventions and the frequency with which policy changes can be made.
The COVID-19 pandemic that emerged from China has caused considerable morbidity and mortality across the globe. Non-pharmaceutical interventions (NPIs), e.g., masking-up in public places, social-distancing, school and border closures, contact-tracing, etc., were crucial in curtailing the burden of the virus during the early stages, while development and use of highly effective vaccines have been useful during the later stages of the pandemic. Despite these non-pharmaceutical and pharmaceutical intervention measures, constraining the pandemic remains challenging in many parts of the world. This is due to several factors that include the emergence of new variants of concern against which existing vaccines are not very efficient, vaccine hesitancy, and low availability of vaccines in some parts of the world. In this study, a mathematical model is developed and used to study the combined impact of pharmaceutical interventions, pharmaceutical interventions, and various variants of concern on the burden of COVID-19 in Cameroon. The model is trained with COVID-19 case and vaccination data from Cameroon. Results of the study indicate that early application of NPIs (specifically masking-up with highly effective masks such as N95 masks) would have prevented the emergence of most of the cases in Cameroon. Additionally, the study shows that herd immunity can be attained if 81% of the population is fully vaccinated, and that this threshold is even higher in the case in which immunity wanes or more transmissible variants of concern are considered. Furthermore, the study shows that striking an appropriate balance between the number of fully vaccinated individuals and the number of individuals who mask-up regularly in public can lead to a drastic decrease in the number of cases in Cameroon.
The evidence for vector control methods aimed at reducing the burden of Aedes-borne diseases largely consists of studies that measure entomological endpoints for a single intervention. In practice, in the effort to control outbreaks, multiple vector control methods are often implemented simultaneously, and it is currently not well understood how different vector control methods interact with each other to reduce disease burden. To address this, we conducted a systematic literature review to obtain estimates of entomological parameters under the impact of eight different vector control methods to calculate transmission coefficients under a Ross-Macdonald formulation. Using the next-generation method, we calculated the reproduction number under the impact of pairs of interventions for a range of coverage levels to determine which combinations resulted in the greatest reduction of transmission. Initial results suggest that as coverage of interventions that increase mortality during early life stages, such as larviciding, increases, interventions that primarily derive their effects from their impacts on vectors during later life stages, such as spatial repellents, exhibit smaller impacts. On the other hand, the impact of interventions that act on overlapping life stages increases as coverage of both interventions increase. Utilizing the next-generation method allows us to effectively investigate potential synergies between pairs of interventions. This method could be extended to exploring synergies between interventions for other infectious diseases.
Infectious disease dynamics in amphibians, which can be transmitted through multiple routes, constitute a complex and interconnected system. The spread of infection varies depending on the population level and age stages of the host species, such as larvae, efts, and adults. Due to seasonal reproductive behaviors and metamorphosis, the population density of each stage fluctuates over time. To study the transmission dynamics of a fungal pathogen, Batrachochytrium Salamandrivorans (Bsal), in North American salamanders across different population densities and environments, we developed a compartmental model using ordinary differential equations. By analyzing model and simulations, we gained insights into strategies for controlling transmission and preventing epidemic outbreaks resulting from different pathogen loads at different temperatures.
In 2022, a resurgence of the cholera outbreak emerged, posing a renewed threat to public health. It can be transmitted through indirect transmission (environment-to-person) by ingesting food or water contaminated with the bacterium Vibrio cholerae. Since climate change is causing shifts in precipitation patterns globally, it can affect the movement of pathogens through stream networks and result in changes in disease dynamics. To investigate the impact of the change, we formulated a multi-patch model for cholera with precipitation level, which affects the stream network. We show the correlation between the basic reproduction and the level of precipitation analytically and numerically. Moreover, we investigated patch-specific optimal treatment strategies.
During the COVID-19 pandemic in past three years, a series of computational and mathematical approaches have been suggested to figure out the epidemic characteristics, which include the effectiveness of social distancing, vaccinations, and the spread itself. In spite of these efforts, high evolution rate of SARS-CoV-2 bears dominant variants of COVID-19 every four to eight months, which leads to failures of improving feasibility of long-term models and understandings. We also witnessed the latest dominant variant Omicron shows a three times higher transmission rate and limits two-dose vaccination against symptomatic infection. We suggest a intergrated mathemathical model, which incorperates three variants of COVID-19 at once, to understand daily pattern from October 2021 to June 2022. It separates subsequent dominant variant occupying twenty percent of the reported cases to GISAID. Indistinguishable patterns are observed in COVID-19 cases from USA, UK, Japan, and in South Korea. We are able to improve the viability of a four months COVID-19 incidence model by dividing the models according to the dominant variant in each period respectively. Based on these, we suggest that consideration of a change of dominant variant of SARS-CoV-2 is necessary in improvement of feasibility of short-term designed stochastic models to a longer-term prediction.
The COVID-19 pandemic has had a significant impact on both overall mortality and COVID-19 deaths worldwide. Estimating excess mortality during the pandemic is a key measure for assessing its direct and indirect effects on public health. To understand the scope of excess mortality during the pandemic in South Korea, we used monthly death and mean temperature data for each region from January 2010 to December 2019 to develop linear models and estimate expected deaths without the pandemic. Our analysis revealed significant regional variation in excess mortality, with some areas experiencing higher rates of excess deaths not attributed to COVID-19. These findings underscore the need for targeted interventions and public health measures to address the indirect effects of the pandemic on mortality, particularly in areas with higher excess mortality. By understanding the patterns of excess mortality and the factors that contribute to regional variation, we can develop more effective strategies to mitigate the impact of the pandemic and protect vulnerable populations.
Bone remodeling is a very complex and fine-tuned process, which is necessary to ensure a healthy bone structure. If this process gets out of balance – e.g., because of hormonal disbalance or the impact of bone metastases – pathologies like osteoporosis can appear. In this contribution we introduce a novel computational approach to investigate this balance by connecting the bone remodeling process with its microenvironment. Our goal is to better understand the well-balanced and complex dynamic of the subprocesses involved in healthy bone remodeling.
We implement a 3D stochastic cellular automaton (SCA), where voxels interact only with their nearest neighbors in a scaffold representing bone tissue. At each time point, each voxel can take one of four different states that stand for the different phases of bone remodeling: formation, quiescent bone, resorption, and environment. To create a compact representation of the frequency-dependent interaction of those voxel states we make use of methods borrowed from evolutionary game theory for the update rule of the cellular automaton [1]. This representation encodes knowledge about the mutual impact the different actors of bone remodeling (osteocytes, osteoclasts and osteoblasts) have on each other. Each parameter in the model has therefore a direct connection to the biological processes. First, we set up simulations of the model with either only resorption or only formation. This choice reduced the model complexity and allowed us to determine parameter spaces for a self-regulating behavior for each of them. The self-regulating behavior is defined by resorption or formation starting and ending without further parameter tuning. Parameters outside the range of self-regulation will lead to either osteolytic lesions (resorption) or heterotopic ossification (formation). Further analyses supported the approach of a spatial model with a small neighborhood to simulate the local phenomena observed in bone remodeling. Next, we coupled the two processes of resorption and formation. In the limit of separation of time scales, our model showed that self-regulating resorption followed by self-regulating formation reproduces the physiological bone remodeling behavior. Further analysis will create a more fluid coupling of the two processes while involving more parameters.
The model has the potential to use the role of the microenvironment to evaluate the impact of additional factors, such as drugs or bone metastases. We are planning on using experimental in vivo data from a breast cancer bone metastasis mouse model [2], which includes spatial and temporal dynamic of early osteolytic lesions, to fit additional parameters. Hopefully, these findings will add to the discussion, how pathological behavior might be controlled, if not even reversed.
[1] M. D. Ryser and K.A. Murgas, Bone remodeling as a spatial evolutionary game, Journal of Theoretical Biology, 2017
[2] S. A. E. Young, A.-D. Heller et al., From breast cancer cell homing to the onset of early bone metastasis: dynamic bone (re)modeling as a driver of osteolytic disease, bioRxiv preprint
Mean-field models of protein translocation in mammalian cell metabolism in response to insulin have previously been used to identify dominant processes at the macroscopic scale (J. Biol. Chem., 289(25): 17280-17298). These mean field models do not take the stochasticity and variance of the data fully into account, however. These models also do not provide explanatory mechanisms for the response to the insulin signal. We have developed a candidate stochastic queuing network model that may provide further insight into mechanisms at the molecular scale for glucose transporter translocation in insulin regulated metabolism.
To test the efficacy of the model as an explanation of the biological mechanisms, an assessment of the ability of the model to represent all the different observations needs to be quantified. For each particular experimental protocol the data set consists of small numbers of repeated samples at discrete time points of the system under that experimental condition. The stochastic model then aims to describe all the different time evolving distributions corresponding to the different experimental protocols.
Not only do the distributions of the data and model need to be compared at each time point in the data set for each protocol, but also a comparison needs to be made across time as each of the distributions evolve. Additionally, the correspondence of the stochastic model and observations across the different experimental protocols needs to be quantified. In systems where data is sparse, robustness can be given to inference when independent data sets from multiple sources are combined, given that the model parameters constrained by the different protocols overlap.
In this investigation, different distance measures and comparators of evolving distributions are explored for the candidate model of glucose transporter translocation with a view to building a practical algorithm for inference of stochastic models with multiple stochastic data sets from different experimental protocols. The efficacy and implications of different approaches and for the candidate model is discussed.
Wet age-related macular degeneration (AMD) is a blinding chronic eye disease commonly treated with monthly intravitreal injections. Drug delivery systems (DDS) aim to reduce injection frequency. Here, we developed mathematical models of drug release from bi-layered prototype chitosan-polycaprolactone (PCL) DDS to help optimize their design and improve wet AMD treatments. Fick’s second law is used to model the unsteady-state drug release from DDS into phosphate buffer saline. For drug-loaded chitosan-PCL microspheres, we solved the diffusion equation numerically using finite differences in MATLAB, and finite elements in COMSOL. We then use COMSOL for modeling a more complicated geometry consisting of a chitosan-PCL cylindrical device with a hollow core for drug loading. Furthermore, we use ordinary least squares objective functions in both software to estimate relevant parameters from the DDS using experimental data. Our MATLAB and COMSOL models accurately simulated the cumulative drug release behavior from the microspheres for 160 days compared to in vitro experimental data. For the cylindrical device, we observed large deviations in the initial 50 days, with more accurate predictions after that, implying other drug-release mechanisms, like erosion, need to be considered for the initial phase. The models can help optimize the design of bi-layered DDS to improve wet AMD treatments and provide insights into the mechanisms involved in the drug release from these DDS.
Transcriptional repression can occur via various mechanisms, such as blocking, sequestration and displacement. For instance, the repressors can hold the activators to prevent binding with DNA or can bind to the DNA-bound activators to block their transcriptional activity. Although the transcription can be completely suppressed with a single mechanism, multiple repression mechanisms are used together to inhibit transcriptional activators in many systems, such as circadian clocks and NF-κB oscillators. This raises the question of what advantages arise if seemingly redundant repression mechanisms are combined. Here, by deriving equations describing the multiple repression mechanisms, we find that their combination can synergistically generate a sharply ultrasensitive transcription response and thus strong oscillations. This rationalizes why the multiple repression mechanisms are used together in various biological oscillators. The critical role of such combined transcriptional repression for strong oscillations is further supported by our analysis of formerly identified mutations disrupting the transcriptional repression of the mammalian circadian clock. The hitherto unrecognized source of the ultrasensitivity, the combined transcriptional repression, can lead to robust synthetic oscillators with a previously unachievable simple design.
Drug discovery is a highly complicated and time-consuming process. One of the main challenges in drug development is predicting whether a drug-like molecule will interact with a specific target protein. This prediction is crucial in expediting the validation and discovery of targets, and it enables biochemists and pharmacists to accelerate the drug development process. In recent studies of biomolecular sciences, the application of algebraic graph-based models to accurately represent molecular complexes and predict drug-target binding affinity has generated significant interest among researchers. Here, we present algebraic graph-based molecular representations to form data-driven scoring functions (SF) named AGL-EAT-Score featuring extended atom types to capture wide-range interactions between the target and drug candidate. Our model applies multiscale weighted colored subgraphs for the protein-ligand complex where the graph coloring is based on SYBYL atom-type and ECIF atom-type interactions. Furthermore, combined with the gradient-boosting decision tree (GBDT) machine-learning algorithm, our newly developed SF has outperformed numerous state-of-the-art models in PDBbind benchmarks for binding affinity scoring power, and the D3R dataset, a worldwide grand challenge in drug design.
Understanding altered metabolism in different conditions requires consideration of various connections across multiple scales. This project aims to understand the metabolic relationship between Colorectal Cancer Cells and Cancer-Associated Fibroblast (CAF). Firstly, an experimental workflow is designed to measure the effect of CAF presence on CRC metabolism. The Flux Balance Analysis (FBA) model is created using growth and metabolomic data. Next, 3-D multiscale agent-based model (ABM) is built to scale from a single cell level to dozens of organoids. We integrated the metabolic model as an FBA model to be employed as a chemical network in each agent. The multiscale model provides the spatial information, which is local substrate availabilities and cellular pressures, to be used as input to FBA. The metabolic model yields a biomass creation rate used as cellular volume growth in agents. Individual agents proliferate with adequate cellular volume and exchange rates for essential chemicals. The distribution of important metabolites in the 3-D domain is calculated by 3-D reaction-diffusion equations. However, the whole computational framework is expensive; therefore, we enhanced our framework with a surrogate model and multiple domains. The metabolic model portion of the simulation is speeded up with a deep neural network (DNN) which is trained by high throughput pre-run FBA model screens. Other acceleration is gained by coarsening the microenvironment domain, which does not contain cells. Multiscale simulations have matched with experimental growth rates. Overall, we combine multiple scales from the molecular level to the 3-D experimental well-containing hundreds of thousands of cells. High-throughput simulations with multiscale knockdowns will help us understand the altered metabolism and discover important targets to diminish this metabolic relationship.
Age-related macular degeneration (AMD) is an irreversible disease caused by macular deterioration and responsible for vision loss. AMD is caused by the growth of abnormal leaky blood vessels due to the high presence of vascular endothelial growth factor (VEGF) in the macular region of the eye. Anti-VEGF drugs have been proven most stable medication in AMD treatment that inhibits the action of vascular endothelial growth factor in the macula. One of the most suggested anti-VEGF drugs used in AMD treatment is Bevacizumab using intravitreal injection. In our study, we developed a 3D spherical region of vitreous for the human and rabbit eye to computationally simulate the pharmacokinetic effect of the intravitreally injected drug Bevacizumab. The model is simulated in COMSOL under time-dependent conditions to observe the spatial drug distribution and calculate the concentration profile in the vitreous and near macula regions. The vitreous is treated as a Darcy porous medium, and the drug transport through the porous medium is solved using mass transport physics coupled with Darcy’s law, including the convection-diffusion effect. The model includes the drug elimination route both anteriorly and posteriorly. Both models are validated against the experimental pharmacokinetic model data using the drug Bevacizumab, and our drug concentration-time plots in vitreous for both the human and rabbit eye are in good agreement with the experimental data. The drug concentration near the macula is also explained with experimental validation.
Population structure has been known to substantially affect evolutionary dynamics. Networks that promote the spreading of fitter mutants are called amplifiers of natural selection, and those that suppress the spreading of fitter mutants are called suppressors. Research in the past two decades has found various families of amplifiers while suppressors still remain somewhat elusive. It has also been discovered that most networks are amplifiers under the birth-death updating combined with uniform initialization, which is a standard condition assumed widely in the literature. In the present study, we extend the birth-death processes to temporal (i.e.,time-varying) networks. For the sake of tractability, we restrict ourselves to switching temporal networks, in which the network structure alternates between two static networks at constant time intervals. We show that, in a majority of cases, switching networks are less amplifying than both of the two static networks constituting the switching networks. Furthermore, most small switching networks are suppressors, which contrasts to the case of static networks.
Corneal injuries following chemical exposure differ in severity and reversibility. Various in vivo, ex vivo, and in vitro experimental methods attempt to predict whether exposure will lead to severe (corrosive), moderate, mild, or no irritancy but differ in their ability to prognosticate human-relevant eye irritation outcome. A detailed computational model of corneal injury at the multi-cellular level (depicting individual cells and biochemical processes in detail) which could predict these adverse outcomes would enable limitless virtual experiments. To improve the spatial and dynamic understanding of corneal chemical hazard, we built a multicellular agent-based model in the CompuCell3D modeling environment that aims to recapitulate complex cell behaviors underlying homeostasis and wound healing of the stratified epithelial layer and the stroma. The model represents a two-dimensional sagittal section of the limbal area with stem and transit-amplifying cells and a stratified epithelium layer keeping the same structure seen in its biological archetype, with a bilayer of superficial cells, two to three layers of wing cells, a single layer of basal cells attached to the basement membrane, and immune cells, bounded by virtual spaces to represent the tear layer and Bowman's membrane. Beneath this epithelial membrane lies an area representing the stroma with keratinocyte cells. Homeostasis in the epithelial layer implements signal information (cytokines, growth factors) and other factors can be added to more completely simulate the emergent wound-healing behavior where tear composition changes after injury, having higher levels of EGF (proliferation and migration), TGF-α (mitogen), HGF (proliferation and migration, promotes wound healing), KGF (proliferation), and IGF (proliferation), in the regulation of composite cellular behavior and multicellular interactions on proliferation and cell migration to the wounded site. These changes in the microenvironment activate quiescent limbal stem cells to proliferate and differentiate into transient-amplifying cells, which also proliferate and consequently differentiate into the other cell types present in the stratified epithelium layer. This mechanism is enough to heal mild and moderate wounds that avoid damaging the basement membrane. In cases of severe injury, other systems, including vascular and myeloid, participate in the repair of the Bowman's membrane and the stroma. This prototype virtual corneal model aims to define a more mechanistic human-relevant classification scheme, predict the time of recovery from each of those injuries, and offer potential explanations for the corneal anomalies (erosions and corneal ulcers) after severe damage and simulated responses to bioactivity data from various in vitro models of corneal toxicity. This will help toxicologists better understand critical events in cornea-chemical exposures as well as predict human-relevant adverse outcomes. Disclaimer: this abstract does not necessarily reflect USEPA policy.
Organisms exhibit a dazzling array of symmetries, from the rotational symmetries of flowers to the fractal symmetries of trees and even bilateral symmetries in humans. Symmetry is fundamental and is often a predictor of survivability, fecundity, and evolvability. Although it is intuitively clear that symmetry exists in nature, the symmetries are typically imperfect, making it difficult to apply mathematical tools that were built to understand idealized versions of symmetry. In 2021, Gandhi et al. proposed a real-valued operator that can quantify approximate symmetries by evaluating how much an object changes under a transformation. When one parametrizes the transformation and considers the operator’s graph on the parameter space, the symmetries of the object appear as local minima. I consider the rotational symmetries of a Turing pattern, showing that if we treat minima and maxima of the graph as stable and unstable equilibria (respectively), the changes in extrema are qualitatively similar to changes in equilibria that we observe in classical local bifurcations. Studying relevant properties of the operator may allow us to apply the tools of bifurcation theory to understand how approximate symmetries form in development.
Biofilms are complex heterogeneous substances that can be viewed from the perspective of soft matter physics and continuum mechanics. Biofilm structure can be modeled as a multiphase system where each component has its own rheological characteristics. From the biophysics point of view, the biofilm components create a gel-mix consisting of a polymeric network (polysaccharides) and fluid solvent. The biological and mechanical interactions between these components govern biofilm physics and its spatial variation. We developed a mathematical model to describe the spatiotemporal organization of the biofilm components as a multiphase system where each volume in space is fractionally occupied by the polymeric network and the fluid solvent. The polymeric network is modeled as a viscoelastic fluid that induces viscoelastic stresses due to the rheological behavior of polysaccharides. This viscoelastic stress is a function of the biofilm viscoelastic properties, which are estimated using a Markov Chain Monte Carlo method based on experimental data. The fluid solvent is modeled as a Newtonian fluid, creating viscous stresses within the computational domain. The dynamics of the phases are governed by the conservation of mass and momentum. Each phase moves with its own velocity, introducing a drag force between the phases that is proportional to the velocity difference between the phases. The motion and interaction of the gel-mix components are formulated as a set of equations in an incompressible Navier-Stokes form. These equations are discretized in integral form for infinitesimal control volumes on a two-dimensional staggered grid. This model helps us understand the motion of the biofilm components and can help future researches elucidate the dynamics of polymeric network that forms the backbone of the biofilm.
Chronic kidney disease (CKD) is a family of kidney diseases with various root causes that lead to eventual kidney failure and are characterized by dysfunction of the glomeruli, the functional subunits of the kidneys where blood is filtered. A glomerulus includes the glomerular filtration barrier (GFB) made of the endothelial layer, basement membrane, podocyte epithelial layer, and glycocalyx. The deterioration of the filtration barrier means that the kidney cannot effectively filter solutes from the capillaries, such as proteins, excess water, and other waste products. The functionality of the GFB is measured by the glomerular filtration rate or the rate at which fluid from the capillaries in the kidney is filtered to be excreted. Assessing glomerular dysfunction during CKD requires quantifying the effect of damage in the anatomical ultrastructure of GFB and the unwanted transport of protein through the GFB. Though various methods of assessing glomerular dysfunction exist, current computational models often neglect the glycocalyx as well as the effect of solute and GFB charge. We use open-source software FEBio (Finite Elements for Biomechanics) to simulate fluid transport in different layers of the GFB. FEBio applies continuum biphasic (fluid dynamics/solid biomechanics) theory to describe viscous fluid interactions with porous-hydrated biological tissues. The biphasic fluid-solid interactions (BFSI) solver in FEBio is used to model structures of the glycocalyx, glomerular basement membrane, porous medium, and fluid-solid interactions through the intricate small channels that form the fenestrated endothelial layer and the GBM. Transport equations describe the movement of fluids and solutes from the blood vessel lumen through the GFB. The anatomical ultrastructural parameters for the proposed model were estimated from high-resolution electron microscopy of the glomerular capillary wall. With the information gathered from the electron microscopy images, a “subunit” consisting of the averaged parameterized features of the filter was used to simulate GFB. In addition, ultrastructural parameters were used to design the 3D fluid domain for the simulation using MATLAB and GIBBON, a dedicated biomechanics add-on. The volumetric domain was exported to FEBio, where material properties, boundary conditions, and an analysis step were included for the model. The conditions of the simulation were analogous to the physiological conditions of the in vivo environment. Our simulations showed the flux of solutes (e.g., albumin, glucose, signaling molecules) through the GFB, which can be used to find the glomerular filtration rate (GFR). We intend to simulate the dynamic effects of biomolecular reactions on kidney ultrastructure as it relates to CKD. We use the model to analyze important dynamic phenomena during disease progression, including the widening of the filtration slit, thickening of the glomerular basement membrane, and detachment of the podocyte food processes. By recreating the human anatomy in a computational platform and applying the correct transport phenomena in each tissue layer, the physiological effects on the transport of solutes and glomerular filtration rate can be determined. Understanding the glomeruli’s fluid transport and chemical and physical interactions is critical to provide insights into human development, disease progression, and wound healing possibilities.
As an animal searching for prey performs a random walk, processes in the animal’s nervous system make decisions that produce a distribution of step lengths. I will describe in detail an Agent-Based Model for how an animal’s nervous system might make these decisions. The “agents,” that I call “Simple Abstract Neurons,” are NON-deterministic generalizations of well-known digital logic gates.
I will use a simple version of the model, with carefully-chosen, made-up parameters to illustrate central concepts. I will give a detailed example of how the model parameters can be adjusted to fit a set of empirical data; in this case, from a diving marine predator: an individual blue shark. The SAN model fits the shark data much more closely than the “best fit” of a theoretical, analytical model.
Theoretical studies suggest that when prey is plentiful, an exponential distribution of step lengths is effective. Otherwise, a power-law distribution is optimum. Animals follow such distributions only to the degree that evolutionary pressures may have resulted in an approximation that gives an acceptable balance of cost vs. benefit. But theory provides no insight as to how an animal might produce the observed behavior. The SAN model suggests some answers and makes testable predictions. For example, the model easily and naturally explains how an animal can follow an approximate power-law distribution while avoiding the implied infinities at very short and very long step lengths.
Because the underlying processes are stochastic, any set of empirical data is a sample from a range of possible sets. Model-generated “synthetic data” can be used to characterize this range. The model can also be used to perform very insightful “What if?” experiments.
SANs are compatible with digital logic; the number of possible pure and hybrid networks is huge. The potential is very large for the SAN model to be adapted to other animal behaviors besides random walks.
Dengue is one of the most common mosquito-borne diseases in the world, and a person can get infected by one of the four serotypes of the virus named DENV-1, DENV-2, DENV-3, and DENV-4. After infection with one of these serotypes, an individual will maintain permanent immunity to that serotype, and partial immunity to the other three serotypes. Therefore, there is a risk of getting infected by this virus a maximum of four times, and the symptoms may vary from mild fever to high fever, bleeding, enlarged liver, and severe shock, and sometimes these symptoms may lead to death. It is obvious that the increase in the number of infected individuals makes a negative impact on a country’s economy. Hence, the use of different control measures such as mosquito repellents and the introduction of a vaccine against the virus is important in controlling the spread of the virus. In this study, I am presenting a methodology on how to estimate the optimal rate of vaccinations based on the QDENGA dengue vaccine and the optimal rate of control measures to reduce the number of new and severe dengue cases while minimizing the overall cost. In addition, this vaccine claims high protection against symptomatic disease and waning protection over time for some DENV serotypes. However, the extent to which protection against disease conditional on infection is unknown. I consider different scenarios subject to the possible combinations of vaccine protection and control measures to investigate the most effective parameter values to control the transmission of the virus. Disease forecasts including the number of newly infected individuals in each serotype, the optimal rate of control measure, and vaccinations for a period of ten years are performed with the help of computer software.
We introduce Catalyst.jl, a flexible and feature-filled Julia library for modeling and high performance simulation of chemical reaction networks (CRNs). Catalyst acts as both a domain-specific language and an intermediate representation for symbolically encoding CRN models as Julia-native objects. This enables a pipeline of symbolically specifying, analyzing, and modifying reaction networks; converting Catalyst models to symbolic representations of concrete mathematical models; and generating compiled code for use in numerical solvers. Currently, Catalyst supports conversion to symbolic discrete stochastic chemical kinetics (jump process), chemical Langevin (stochastic differential equation), and mass-action reaction rate equation (ordinary differential equation) models. Leveraging ModelingToolkit.jl and Symbolics.jl, Catalyst models can be analyzed, simplified, and compiled into optimized representations for use in a broad variety of numerical solvers. The performance of the numerical solvers Catalyst targets is illustrated across a variety of reaction networks by benchmarking stochastic simulation algorithm and ODE solver performance. These benchmarks demonstrate significant performance improvements compared to several popular reaction network simulators. Finally, Catalyst combines with a range of packages within the Julia package ecosystem, enabling functions such as steady state finding, bifurcation analysis, parameter fitting, and much more.
Beta cells are cells in the pancreas that produce and release insulin in response to blood glucose levels. Interactions between beta cells within their local network of an islet is important for the regulation of insulin secretion and to enhance the glucose stimulated response. Beta cells are coupled through gap junctions and generate synchronous threshold-based oscillations of their membrane potential. Dysfunction of coupling has been associated with diabetes. Experiments have suggested individual beta cells can control synchronization. We have previously shown in specific conditions a 'switch' cell can serve this purpose. However, the cellular and network conditions are not fully understood. To test a minimal model representation of this behavior we use a mathematical model of bursting in two triplet configurations, chain and triangle. Biological heterogeneity is introduced by varying the gap junctional coupling and the rate of calcium extrusion parameters for each of the cells, which permits varying types of frequencies. We measure the amplitude of a patched steady cell, and we investigate how the bursting of a high frequency cell and coupling can lead to change of the behavior of the patched steady cell. To demonstrate a switch cell exists, we effect the second intermediate frequency cell by (a) silencing it setting the voltage to rest or (b) ablating it disconnecting this cell from other cells, and observe under what conditions there is a loss of activity. We have found the range of coupling strength and calcium extrusion parameters that support switch cell behavior in the simplified system.
We will present methods to uncover the unknown dynamics and features of several biological systems via deep neural network (DNN). We will show how DNNs can be used as approximations to flow maps of the true underlying biological system utilizing residual networks. Further, we will demonstrate its extension to systems with only partially observed data and systems with uncertain parameters. Once an accurate DNN model is constructed, it can be used as a predictive model for the unknown system, allowing us to conduct further system analysis. We will further explore its extension to the inverse problem, or recovering parameters on the system given the available data.
We apply a reduction method to a published uterine smooth muscle cell (uSMC) by Tong, et al. 2011, using both a representative ion channel and steady-state approximation approach. Although an extensive catalogue of potassium channels are known to reside in the uSMC, we hypothesise not all are functionally relevant to reproduce the data given. Further, the Tong model incorporates a vast range of time scales for the Hodgkin-Huxley type activation and inactivataion variables ranging over six orders of magnitude. We demonstrate effective use of a reduced suite of potassium channels and deployment of steady-state approximations that not only reproduces the same data set as the Tong model but additional data a later expanded Tong model with even more potassium channels. Moreover, our reduced model increases computational performance by 200%.
The circadian clock is an internal timer that coordinates the daily rhythms of behavior and physiology, including sleep and hormone secretion. Accurately tracking the state of the circadian clock, or circadian phase, holds immense potential for precision medicine. Wearable devices present an opportunity to estimate the circadian phase in the real world, as they can non-invasively monitor various physiological outputs influenced by the circadian clock. However, accurately estimating circadian phase from wearable data remains challenging, primarily due to the lack of methods that integrate minute-by-minute wearable data with prior knowledge of the circadian phase. To address this issue, we propose a framework that integrates multi-time scale physiological data to estimate the circadian phase, along with an efficient implementation algorithm based on Bayesian inference and a new state space estimation method called the level set Kalman filter. Our numerical experiments indicate that our approach outperforms previous methods for circadian phase estimation consistently. Furthermore, our method enables us to examine the contribution of noise from different sources to the estimation, which was not feasible with prior methods. We found that internal noise unrelated to external stimuli is a crucial factor in determining estimation results. Lastly, we developed a user-friendly computational package and applied it to real-world data to demonstrate the potential value of our approach. Our results provide a foundation for systematically understanding the real-world dynamics of the circadian clock.
Smelling is key for many cognitive and behavioral tasks and is particularly unique having two modes: through the nasal cavity from the front (sniffing) or rear (eating), i.e., orthonasal and retronasal, respectively. Little is known about the differences in how olfactory bulb (OB) cells process odor information between these two modes (ortho/retro). Based on multi-electrode array recordings in rat OB, we find significant differences between ortho and retro spiking statistics – the mode (ortho/retro) is encoded. Using GABA_A agonists and antagonists, we find intermediate levels of inhibition give the best average decoding accuracy of ortho vs retro odors. Our computational models show how inhibition effects decoding accuracy.
Drug sensitive cancer cells often acquire drug resistance, resulting in cancer recurrence despite an initial reduction in tumor size. Two common mechanisms for acquiring drug resistance are point mutation and gene amplification. We propose stochastic multi-type branching process models for each of these mechanisms. Using these models, we derive tumor extinction probabilities and deterministic estimates for the tumor recurrence time, that is, the time when an initially drug sensitive tumor surpasses its original size after developing resistance. For each model, we prove a law of large numbers result regarding the convergence of the stochastic recurrence time to its mean. Additionally, we prove sufficient and necessary conditions for a tumor to escape extinction under the gene amplification model, discuss behavior under biologically relevant parameters, and compare the recurrence time and tumor composition in the mutation and amplification models both analytically and using simulations.
Tumors develop in a complex physical, biochemical, and cellular milieu, referred to as the tumor microenvironment. Of special interest is the set of immune cells that reciprocally interact with the tumor, the tumor-immune microenvironment (TIME). The diversity of cell types and cell-cell interactions in the TIME has led researchers to apply concepts from ecology to describe the dynamics. However, while tumor cells are known to induce immune cells to switch from anti-tumor to pro-tumor phenotypes, this type of ecological interaction has been largely overlooked. To address this gap in cancer modeling, we develop a minimal, ecological model of the TIME with immune cell conversion, to highlight this important interaction and explore its consequences. A key finding is that immune conversion increases the range of parameters supporting a co-existence phase in which the immune system and the tumor reach a stalemate. Our results suggest that further investigation of the consequences of immune cell conversion, using detailed, data-driven models, will be critical for greater understanding of TIME dynamics.
Neural stem cells (NSCs) have been shown to be a promising treatment for various brain pathologies due to their ability to migrate directly to the target site, repair damaged tissue, and deliver therapeutic agents. However, the efficacy of such treatments relies on the number and timing of viable cells that are able to reach the injury site. These factors are greatly influenced by different injection strategies as well the complex cytokine dynamics within the brain. For instance, intracranial injections are highly invasive but can be administered next to the site, while intranasal injections can be administered multiple times but require the cells to travel from the olfactory bulb. Furthermore, NSC’s sensitivity to chemoattractants within the brain can alter the path taken and allow for more robust migration. To understand and test the mechanisms that govern NSC migration in a cost effective manner, we build a probabilistic model which accounts for crossing white matter tracts and chemotaxis within the settings of naive rat/mouse brain and TBI. We then use the model to predict the migratory paths in response to different injection strategies/timings as well as explore possible combination therapies to increase NSC arrival rates.
While evolutionary approaches to medicine hold great potential, measuring evolution is difficult due to experimental constraints and the dynamic nature of biology. It is impossible to continuously observe the evolution of cancer, and obtaining multiple longitudinal samples over time is rare. Advancements in single-cell DNA sequencing have allowed for new evolutionary approaches to studying somatic clonal expansion, which are likely to improve mechanistic understanding of cancer and our ability to effectively prognosticate patients. We present coalescent methods to estimate the growth rate of clones from reconstructed evolutionary trees, eliminating the need for complex simulations. We apply our methods to four recently published single-cell whole genome sequencing datasets, estimating the growth rate of clonal expansions in blood, and validating these estimates with longitudinal data. We show that our estimates lead to new insights on evolutionary parameters, which have implications for early detection of high-risk clones. For example, compared to clones with a single or unknown driver mutation, clones with multiple drivers have increased growth rates (median 0.94 vs. 0.25 per cell per year; p = 1.6 x10^-6). Additionally, patients diagnosed with Myeloproliferative Neoplasm (MPN), a group of malignant conditions characterized by overproduction of blood cells, were found to harbor more aggressively expanding clones (median 0.55 vs. 0.23 per cell per year; p = 0.029) compared to healthy individuals. Further, stratifying patients with MPN by the growth rate of their fittest clone uncovered that higher growth rates are associated with shorter time from clone initiation to MPN diagnosis (median 13.9 vs. 26.4 months; p = 0.0026). As genomic sequencing technology continues to advance, we demonstrate that clonal growth rates can be accurately estimated and have potential for clinical application. To make our methods widely available, we created cloneRate, a user-friendly R package for researchers to apply to their own datasets.
Lung cancer is a major health problem, with high incidence and mortality rates, due to an absence of effective treatment strategies. In the United States, it is the leading cause of death by cancer, with a 5-year survival of 23% (ACS source). The molecular basis of this disease is complex and heterogeneous, and large inter-patient variability is observed in treatment response, so that it is necessary to consider a mathematical approach to study the processes involved in lung cancer and personalize therapies. In this work, we focus on two deregulated mechanisms in cancer: the immune system and the circadian clock. At the cell level, the circadian clock is a 24-hour biological oscillator that regulates most intracellular processes. It consists of a regulatory network with several intertwine feedback loops that generate sustained oscillations with a period between 20 and 30h. Besides, NLRP3 is a sensor of innate immunity whose role in the immune response has been well studied which was recently identified as an interesting gene altered in lung tumors and predictive of poor prognosis. Previous studies have shown that NLRP3 transcription is regulated indirectly by REV-ERB α (a nuclear receptor of the circadian clock) in macrophages [1][2]. However, the links between NLRP3 and the circadian clock and in particular the impact of the clock on the function of NLRP3 have been very rarely investigated. Our goal is then to characterize the interactions between NLRP3 and the circadian clock that are emerging as major components in the pathophysiology of lung cancer. To this end, we have undertaken a combined experimental and mathematical approach. We have studied the interactions of NLRP3 and the circadian clock in human bronchial epithelial cells (HBEC) which were synchronized by serum shock. Transcriptomics (RNA-Sequencing) and proteomics (Mass Spectrometry) data as well as intracellular localization (nucleus/cytoplasm) were assessed. Clock gene components were defined using the Reactome database (v84). Circadian rhythms were studied using cosine wave fitting and using CMAES for the minimization task. Model learning method was developed to automatically learn the structure of quantitative systems biology models based on ordinary differential equations from multimodal data. Parameter estimation was performed using a modified least-square approach using CMAES for minimization. The analysis of mRNA levels of 67 clock genes revealed a functional clock in HBEC cells with a period of 30h+/-2h . NLRP3 can interact with clock proteins and the data suggest that they could regulate the intracellular localization of NLRP3 to orchestrate its functions. On the other hand, loss of NLRP3 expression may disrupt the circadian regulation necessary for normal lung function. An existing circadian clock model [3] using ordinary differential equations (ODE) was extended by adding equations describing the influence of the clock on NLRP3 transcription and interactions of clock and NLRP3 proteins. As a start, a collection of models were considered that included a single additional reaction as compared to the initial clock model. Datasets used for the fit were: mRNA levels of 7 clock genes, protein level of 7 clock genes and circadian rhythms of nucleus/cytoplasm localization of NLRP3, BMAL1, PER2 and CRY1. A systematic fit of each model was performed which allowed to eliminate unlikely reactions. Models involving more than one additional reaction are being investigated. Such model learning pipeline will help prioritize future experiments to fully determine NLRP3 interactions with the clock and identify potential drug targets to restore NLRP3 functions in NLRP3-altered cancer cells.
[1] Pourcet, B., Zecchin, M., Ferri, L., Beauchamp, J., Sitaula, S., Billon, C., ... & Duez, H. M. (2018). Nuclear receptor subfamily 1 group D member 1 regulates circadian activity of NLRP3 inflammasome to reduce the severity of fulminant hepatitis in mice. Gastroenterology, 154(5), 1449-1464. [2] Wang, S., Lin, Y., Yuan, X., Li, F., Guo, L., & Wu, B. (2018). REV-ERBα integrates colon clock with experimental colitis through regulation of NF-κB/NLRP3 axis. Nature communications, 9(1), 1-12. [3] J. Hesse, J. Martinelli, O. Aboumanify, A. Ballesta, and A. Relogio. A mathematical model of the circadian clock and drug pharmacology to optimize irinotecan administration timing in colorectal cancer. Computational and structural biotechnology journal, 19:5170–5183, 2021.
Chronic myeloid leukemia (CML) is a blood cancer characterized by dysregulated production of maturing myeloid cells driven by the product of the Philadelphia chromosome, the BCR-ABL1 tyrosine kinase. Tyrosine kinase inhibitors (TKI) have proved effective in treating CML but there is still a cohort of patients who do not respond to TKI therapy even in the absence of mutations in the BCR-ABL1 kinase domain that mediate drug resistance. To discover novel strategies to improve TKI therapy in CML, we developed a nonlinear mathematical model of CML hematopoiesis that incorporates feedback control and lineage branching. Cell-cell interactions were constrained using an automated model selection method together with previous observations and new in vivo data from a chimeric BCR-ABL1 transgenic mouse model of CML. The resulting quantitative model captures the dynamics of normal and CML cells at various stages of the disease, exhibits variable responses to TKI treatment, predicts key factors of refractory response to TKI treatment, and predicts potential combination therapy efficacy. Recent experiments reveal that interactions and competition between different cellular compartments and between normal and BCR-ABL1-expressing cells form a threshold that determines whether the malignant cells can expand and cause leukemia. To capture these experimental dynamics, we found it necessary to incorporate additional biological factors through the introduction of new cell types and interactions. We applied an adapted model selection scheme to explore the unknown cell-cell interaction space and find subsets of models consistent with experimental dynamics. We analyzed common motifs across experimentally consistent models and identified interactions as targets for experimental design to further narrow the valid models.
In recent decades, research on the active expression and regulatory effects of chemokines in cancer and immune cells has made the chemokine system an emerging target of immunotherapy. Alteration in chemokine environments is expected during immunotherapy, emphasizing the importance of understanding cell migration in complex chemokine environments. The complex signaling network formed by chemokines and cognate receptors regulates diverse tumor and immune cell activities, including leukocyte recruitment, angiogenesis, tumor growth, proliferation, and metastasis. We built 2D & 3D agent-based models with Compucell3D (a cellular Potts lattice-based model) to simulate the physiological response, especially cell migration, of tumor and immune cells towards complex chemokine settings. The 2D model is used to understand the mechanisms of cell chemotaxis, monomer-dimer equilibrium of certain chemokines, and competition between different pairs of chemokines and cognate receptors. The 3D model simulates and predicts an in vitro transwell experiment where cells have more realistic biomechanics of neighboring cells and tissue-mimic biomaterials. Using the models, we investigated how chemokine concentration, chemotactic force, environment composition, energy term that governs random walk, and membrane properties can influence cell migration. Results from this study will be used to build new agent-based models to simulate in vivo cancer pathology and therapy, considering cells, chemokines, and tissue microenvironments.
Cancer metastasis accounts for many cancer-related deaths worldwide. Metastasis occurs when the primary tumor sheds cells into the blood and/ or lymphatic circulation, thereby becoming circulating tumor cells (CTCs) that transverse through the circulatory system, extravasate the circulation and establish a secondary distant tumor. CTCs transition through the blood system, which has a plethora of supportive and antitumoral immune cells, represents one of the major metastatic hurdles that are not yet fully deciphered. Upon their entry into the blood stream, CTCs interact with platelets which shield them from the recognition by immune cells, including natural killer cells and CD8+ T cells. Platelet binding to CTCs also enhances CTC arrest in the vascular endothelial walls and subsequent extravasation. On the other hand, activated circulating CD8+ T cells, are able to recognise and attack the arrested CTCs prior to their extravasation. Thus, understanding how the dynamic interactions between CTCs, platelets and CD8+ T cells eventually result in secondary metastatic tumor emergence is a key challenge. Here, through a simple mathematical model of ordinary differential equations, we provide our perspective on how CTCs mechanistically evade the CD8+ T cell cytotoxicity. To achieve this, we aim to (a) describe how the intrinsic growth of the primary tumor, and subsequent dissimination of CTCs, link to the secondary establishment of distant tumor cell population; (b) illustrate how the intravascular proliferation of arrested CTCs within the circulation contributes to the possibility of secondary metastasis; (c) describe possible mechanisms underlying the antitumoral activity of CD8+T cells in inhibiting metastatic potential of CTCs; (d) discuss how simple treatment scenarios can be employed to minimize a further spread of CTCs within the circulation, by focusing on the CTCs disseminated from the primary tumor, rather than the secondary metastatic tumor. Last, we also provide a comprehensive mathematical stability analyses to assess different treatment scenarios that can hamper CTCs survival and highlight the significant role of mathematical modeling in clinical oncology.
BACKGROUND: Acute myeloid leukemia (AML) outcomes remain poor, likely due to treatment-resistant leukemic stem cells (LSCs). Evidence suggests resistance to tyrosine kinase inhibitors (TKIs) depends on the bone marrow’s vascular plasticity via the Janus Phenomena. The opposing “faces” of the Janus Phenomena are the initial beneficial cytoreduction of blast cells, followed by revascularization of the BM by endothelial cells, the expansion of LSCs, and then relapse. This is partially mediated through miR-126 over-expression upon which the endothelial cells are dependent to revascularize the BM and to shelter LSCs. To mitigate the Janus Phenomena, miRisten, a miR-126 expression inhibitor, was developed. To understand the interactions of TKI with miRisten, we used experimental data to inform an agent-based model (ABM) recapitulating the Janus Phenomena to explore how different vasculature architectures and drug scheduling can prevent relapse. In addition to the BM structure, we explored common interactors among genes linked to both drugs (AC220, a TKI, and miRisten). METHODS: Using an on-lattice 2D ABM containing three agents (EC cells, LSCs, and blast cells), we modeled 16 different vascular architectures informed from mouse tibias post-TKI treatment. We tested three conditions: TKI, TKI + miRisten, and TKI + miRisten pre-treatment. K-means clustering and PCA were performed to determine relationships among the varying vascular architectures. A protein-protein interaction (PPI) network based on publicly available data for the proteins affected by both drugs was constructed. We used several centrality measurements to determine the nodes in the directed graph that have the largest role in the connectivity of the network. Next we constructed a graph neural network to classify proteins linked to either the targets of AC220 or miRisten with limited or tenuous experimental validation. RESULTS: The optimal dose strategy to prevent relapse was two-weeks of pre-treatment with mRristen before TKI administration. Additionally, vascular architectures spanning the entire domain of the ABM consistently prevented relapse during treatment. Interestingly within the drugs’ PPI, miR-206 (a known tumor suppressor linked to angiogenesis and an indirect regulator of miR-126) was the node with the highest degree centrality. CONCLUSION: To prevent AML relapse under TKI, miRisten should be given for two weeks before starting TKI. Additionally, both the physical structure of the vasculature and the protein-protein interactors contribute to resistance.
Introduction: Predictive mathematical models of glioma growth and therapy response can be informed with quantitative magnetic resonance imaging (MRI) data [1]. However, it is unclear what quantity and quality of longitudinal MRI data are required for accurate model calibration and prediction. To address this uncertainty, we utilize a novel in silico framework that explores the predictability limits of a spatiotemporal reaction-diffusion glioma model by quantifying how different combinations of signal-to-noise ratio (SNR), spatial resolution (SR), and temporal resolution (TR) in initial murine MRI data affect accuracy of parameter calibration and tumor growth prediction. Methods: We have developed a two-species reaction-diffusion glioma model that describes the spatiotemporal evolution of tumor cellularity and vascularity [2]. The initial cellularity and vascularity conditions that inform the model are estimated from quantitative MRI data. Different combinations of SNR, SR, and TR are applied to the initial conditions. TR is varied by changing the quantity and spacing of the MRI data used to inform parameter calibration. We apply the model to a spatially heterogeneous rat tumor in a simulated brain tissue domain [3]. We then solve the model via the finite difference method and calibrate model parameters using the Levenberg-Marquardt algorithm (N = 50 in silico replicates). We report the mean and standard deviation of each model parameter error as well as error in longitudinal tumor volume prediction for each tested combination of SNR, SR, and TR. Results and future directions: Low SR (voxel volume 0.500 mm3), experimentally relevant SR (voxel volume 0.063 mm3), and high SR (voxel volume 0.008 mm3) conditions are evaluated across a range of TR and SNR for all model parameters. The worst TR case uses two calibration timepoints, 96 hours apart; the experimentally relevant case uses three timepoints, 48 hours apart; and finally, the best case uses five timepoints, 24 hours apart. SNR is tested from 5 to 160. At an experimentally relevant SNR of 40, calibrated parameter percent error (PE) in tumor diffusion and proliferation falls by approximately 75% as SR improves. PE falls by approximately 37% as TR improves in this same SNR case. With SR held constant at the 0.063 mm3 voxel volume SR condition, PE decreases by approximately 40% between an SNR of 20 and 40, with little improvement seen past an SNR of 80. Global concordance correlation coefficients and Dice similarity coefficients were relatively consistent across all tested combinations. The next step will be to evaluate the model with different ground truth tumors in order to recommend target combinations of SR, TR, and SNR for a wider range of tumor types and experimental conditions. Funding: CPRIT RR160005, RP220225; NIH R01CA235800, U24CA226110, U01CA174706. References: [1] Hormuth et al., Advanced Drug Delivery Reviews, 187(114367), 2022. [2] Hormuth et al., Cancers, 13(8), 2021. [3] Hormuth et al., Ann. Biomed. Eng., 47(7), 2019.
Recurrent prostate tumors are commonly treated with a total androgen blockade via chemical castration. In turn, cancerous cells have been known to respond with an evolutionary up-regulation of androgen receptors (AR), thus prolonging cell proliferation and delaying apoptosis. Prostate epithelial cancers treated with androgen ablation therapy invariably become castration resistant. However, the primary mechanism remains unknown. Suggested hypotheses include phenotypic plasticity and natural selection. Here we show that castration resistance in prostate cancers treated with androgen ablation arises through natural selection acting on phenotypic plasticity. We found that tumor aggressiveness, measured as growth rate of serum concentration of prostate specific antigen (PSA), correlates positively with the number of treatment cycles. Additionally, we found a signal of increasing tumor aggressiveness with cycle in both on and off-treatment phases. This result argues against the plasticity hypothesis and is consistent with evolution by natural selection. If plasticity were the mechanism, then tumor aggressiveness would not correlate with cycle. This result can help inform clinical management of prostate cancer treated with androgen ablation. Identification of the exact evolutionary mechanism will almost certainly yield insight into more efficacious treatment schedules, and drug combinations, while maintaining patient quality of life and delaying the onset of castration resistance.
An active area of GBM research is identifying how properties of the brain tissue around the tumor as detected by MRI impact tumor growth and treatment response. Such information could be useful in determining surgical margins and optimizing patient specific treatment selection.
To answer this question, we developed a tumor margin propagating algorithm in which the margin growth or shrinking rate can be locally accelerated or slowed according to T1, T1-post, T2, ADC, and FLAIR scan values. We measured how well the growing/shrinking tumor margin starting from a patient scan overlaps with the tumor margin from the subsequent patient scan with these local rate adjustments.
We applied this approach to MRI sequences from 32 patients treated with hypofractionated stereotactic radiotherapy, bevacizumab and pembrolizumab at Moffitt Cancer Center. We found that the tumors tend have affinity for high-FLAIR regions in most cases, and that tumors that are growing tend to have affinity for high-T1 regions, and tumors that are shrinking tend to avoid high-T1 regions. Once our findings from this project are fully developed, they may assist future modeling efforts to predict tumor proliferation and response to immunotherapy in GBM.
Predictive modeling in oncology is a growing field. The calibration of mathematical model parameters based on limited clinical data is critical to reliable predictions per-patient basis. One omnipresent mathematical model is the reaction-diffusion equation, which has been shown to simulate and predict clinical parameters in different cancer types. Here, we focus on analyzing cell-level data routinely obtained from tissue biopsies at diagnosis for most cancers. We analyze the spatial architecture in biopsy tissues stained with multiplex immunofluorescence. We derive the two-point correlation function and the corresponding spatial power spectral distribution. We show that the data-deduced spatial power spectral distribution can fit the spatial power spectrum of the solution of the reaction-diffusion equation, thereby identifying patient-specific tumor growth and invasion rates from a single, routinely collected clinical tissue. This novel approach is essential for model-parameter-inference for tumor infiltration, which may ultimately be used to inform clinical treatments.
The overall purpose of this study is to determine how different assumptions of radiotherapy efficacy affect predictions of tumor cell count using a biology-based mathematical model describing the spatiotemporal evolution of tumor growth and response to radiotherapy. Models seeking to predict patient-specific response have yet to characterize intratumoral heterogeneity in response to radiotherapy. To address this limitation, we acquired quantitative magnetic resonance imaging (MRI) data on four patients with high-grade gliomas at the MD Anderson Cancer Center being treated with fractionated radiotherapy. This longitudinal data was then used to inform a two-species mechanically-coupled reaction diffusion model [1] describing the spatiotemporal change of tumor growth and response to therapy. Tumor cell proliferation rates, tumor diffusion coefficients, and response to radiotherapy (estimated as the surviving fraction following a single radiotherapy session) were calibrated from data up to 1-month post-radiotherapy using the Levenberg-Marquardt approach in MATLAB. With these patient-specific calibrated parameters, our model simulated tumor growth and response assuming treatment efficacy varies homogeneously (globally) or heterogeneously (as a function of vasculature and cell density). We calculated and compared the percent change in tumor cell count three months after initial treatment for surviving fractions of 0.2 to 1 (in increments of 0.05) for four patients. Treatment response as observed at 3-months post-radiotherapy varied greatly (from eradication to residual disease) depending on each assumption on the spatial variations in efficacy. For example, a surviving fraction of 0.6 resulted in complete eradication of the tumor under both homogenous and heterogenous (cell density) assumptions. However, when radiotherapy efficacy was related to vasculature, only an average 55% decrease in tumor cell count was observed. Thus, we have developed an approach to quantify the impact of different assumptions of heterogeneity in response to radiation on percent change in tumor cell count. Future efforts will extend this approach to a larger cohort of patients.
Immunotherapy has dramatically transformed the cancer treatment landscape. Of the variety of types of immunotherapies available, immune checkpoint inhibitors (ICIs) have gained the spotlight. Although ICIs have shown promising results for some patients, the low response rates in many cancers highlight the challenges of using immune checkpoint blockade as an effective treatment. Cytotoxic T lymphocytes (CTLs) execute their cell-killing function via two distinct mechanisms. The first process is fast-acting and perforin/granzyme-mediated, and the second is a slower, Fas ligand (FasL)-driven killing mechanism. There is also evidence suggesting that the preferred killing mechanism by CTLs depends on the antigenicity of tumor cells. To determine the key factors affecting responses to checkpoint blockade therapy, we constructed an ordinary differential equation model describing in vivo tumor-immune dynamics in the presence of active or blocked PD-1/PD- L1 immune checkpoint. Specifically, we analyzed which aspects of the tumor-immune landscape affect the response to ICIs with endpoints of tumor size and composition in the short and long term. By generating a virtual cohort with heterogeneous tumor and immune attributes, we also simulated the therapeutic outcomes of immune checkpoint blockade in a largely diverse population. In this way, we identified key tumor and immune characteristics that are associated with tumor elimination, dormancy and escape. Our analysis sheds light on which fraction of a population potentially responds well to ICIs and ways to enhance therapeutic outcomes with combination therapy.
Introduction: Predictive mathematical models of glioma growth and therapy response can be informed with quantitative magnetic resonance imaging (MRI) data [1]. However, it is unclear what quantity and quality of longitudinal MRI data are required for accurate model calibration and prediction. To address this uncertainty, we utilize a novel in silico framework that explores the predictability limits of a spatiotemporal reaction-diffusion glioma model by quantifying how different combinations of signal-to-noise ratio (SNR), spatial resolution (SR), and temporal resolution (TR) in initial murine MRI data affect accuracy of parameter calibration and tumor growth prediction. Methods: We have developed a two-species reaction-diffusion glioma model that describes the spatiotemporal evolution of tumor cellularity and vascularity [2]. The initial cellularity and vascularity conditions that inform the model are estimated from quantitative MRI data. Different combinations of SNR, SR, and TR are applied to the initial conditions. TR is varied by changing the quantity and spacing of the MRI data used to inform parameter calibration. We apply the model to a spatially heterogeneous rat tumor in a simulated brain tissue domain [3]. We then solve the model via the finite difference method and calibrate model parameters using the Levenberg-Marquardt algorithm (N = 50 in silico replicates). We report the mean and standard deviation of each model parameter error as well as error in longitudinal tumor volume prediction for each tested combination of SNR, SR, and TR. Results and future directions: Low SR (voxel volume 0.500 mm3), experimentally relevant SR (voxel volume 0.063 mm3), and high SR (voxel volume 0.008 mm3) conditions are evaluated across a range of TR and SNR for all model parameters. The worst TR case uses two calibration timepoints, 96 hours apart; the experimentally relevant case uses three timepoints, 48 hours apart; and finally, the best case uses five timepoints, 24 hours apart. SNR is tested from 5 to 160. At an experimentally relevant SNR of 40, calibrated parameter percent error (PE) in tumor diffusion and proliferation falls by approximately 75% as SR improves. PE falls by approximately 37% as TR improves in this same SNR case. With SR held constant at the 0.063 mm3 voxel volume SR condition, PE decreases by approximately 40% between an SNR of 20 and 40, with little improvement seen past an SNR of 80. Global concordance correlation coefficients and Dice similarity coefficients were relatively consistent across all tested combinations. The next step will be to evaluate the model with different ground truth tumors in order to recommend target combinations of SR, TR, and SNR for a wider range of tumor types and experimental conditions. Funding: CPRIT RR160005, RP220225; NIH R01CA235800, U24CA226110, U01CA174706. References: [1] Hormuth et al., Advanced Drug Delivery Reviews, 187(114367), 2022. [2] Hormuth et al., Cancers, 13(8), 2021. [3] Hormuth et al., Ann. Biomed. Eng., 47(7), 2019.
In this poster, I will summarize the content of my PhD thesis. The main concern of my studies has been mathematical modelling of in vitro cancer cell migration. Along with themodelling, an array of Bayesian statistical inference algorithms for key parameters in the models are presented. The guiding principle behind my research interest is that solid models derived from physical principles can aid in the understanding of how cancer cells interact with one another. The subsequent clinical applications of this research can for example be profiling of cells sampled from a specific patient, aiding the physician in choice of clinical intervention. My model paradigm of choice are agent-based models, where every single cell in the sample is given consideration as an agent. The fundamental building block is a set of stochasticdifferential equations (SDE:s) describing the current location of all cells. We also incorporate cell proliferation into our model, every cell divides or die according to a non-homogeneous Poisson process depending the state of the population.
SMB2023 Follow
Thursday, July 20 at 6:00pm during the "PS02" time block. Room assignment: Performance Hall.