Cardiac allograft vasculopathy (CAV) is a coronary artery disease affecting 50% of heart transplant recipients, and it is the major cause of chronic graft rejection. CAV is driven by the interplay of immunological and non-immunological factors with the infiltration of macrophages as one of the main pathological triggers, setting off a cascade of events promoting endothelial damage and vascular cell dysfunction. Since etiology and evolution of the pathology are still largely unknown, disease management remains challenging and re-transplantation is today the only long-term solution to CAV. A deep understanding of the pathology mechanobiology is fundamental to improve prevention, diagnosis, and treatment of CAV. So far, in vivo models, mostly mouse-based, have been widely used to study CAV, but they are resource-consuming, pose many ethical issues, and allow limited time points of investigation during experimental follow-up. Recently, agent-based models (ABMs) proved to be valid computational tools for capturing and deciphering processes at cell/tissue level, augmenting cost-effectively in vivo lab-based experiments, i.e., guaranteeing richness in observation time points while maintaining low resource consumptions. We hypothesize that integrating ABMs with lab-based experiments can aid classic pre-clinical research by overcoming its limitations. Accordingly, we present a bidimensional ABM of CAV in a mouse-like coronary artery cross-section, simulating the arterial wall response to two distinct stimuli: inflammation and hemodynamic disturbances, the latter in terms of low wall-shear stress (WSS), which together trigger macrophage response activation and exacerbate vascular cell activities. In addition, we performed an extensive analysis to investigate the ABM working mechanisms and gain insight on the driving parameters and the stimuli influences. The ABM replicates with high fidelity a 4-week CAV initiation and progression, well highlighting lumen area decreasing due to progressive intimal thickening in regions exposed to high inflammation and low WSS. The sensitivity analysis remarked that the inflammatory-related events, rather than the WSS, predominantly drive CAV, corroborating the inflammatory nature of the vasculopathy. This proof-of-concept model offers to the scientific community an agile computational platform to deepen CAV understanding and to support the in vivo analysis of CAV in a cost-effective fashion.

The mechanisms employed by blood vessels in order to adapt to their hemodynamic environment are important for our general understanding of disease and development. Changes in arterial geometry are generally induced by two effects: vasodilation and/or constriction; and growth and remodeling (“G&R”). The first can occur over short periods of a few minutes, while the second usually occurs over timescales of weeks or months. The free radical Nitric oxide (NO) is one of the few biological signaling molecules that is gaseous. When smooth muscle cells internalize NO, they lengthen and ultimately induce a relaxation of the artery. In addition, Platelet-Derived Growth Factor (PDGF) is a growth factor released by smooth muscle cells and platelets that regulates cell growth and division. In this talk, we present a single-layered, axisymmetric hyperelastic model for a deforming, growing artery in which the opening angle is regulated by NO and growth is induced by PDGF. Our model describes vasodilation and G&R in a long cylindrical artery regulated by a steady-state Poiseuille flow. The transport of NO released by the endothelium is governed by a diffusion equation with a shear-stress dependent flux boundary condition. Arterial opening angle is assumed to be a Hill function of the wall-averaged NO concentration. We find that both growth and NO help to regulate shear stress with respect to the flow rate, but regulation through growth occurs only at large times. In contrast, regulation through NO is immediate but can only occur as long as the opening angle is able to continually decrease as a function of flow rate. Our model is calibrated using experimental data from ligated, control, and anastomosed carotid arteries of adult and weanling rabbits. Our results generate shear stress/flow rate and lumen radius/flow rate curves that agree with experimental data from control and NO-inhibited rabbit carotid arteries.

Cardiovascular disease (CVD) is presently one of the leading causes of death, with an estimated 24.1 million people expected to be affected by 2025. Therefore, the establishment of the health care industry's objective is to gather a vast amount of data on cardiovascular disease and utilize Deep Learning (DL) algorithms to analyze the information to assist doctors in early detection and identification of potential risk factors for CVD. DL algorithms can help to discover potential patterns of diseases and symptoms based on this structured and unstructured case information. In epidemiology, this is the first prospective study on cardiovascular disease in the community free movement population, and the related risk factors can be recognized. The prediction method of cardiovascular disease based on LSTM is proposed, and the connection between LSTM and unit state is tried to ensure the correct data acquisition during operation, and the prediction method based on LSTM is realized. The original medical data of 4434 participants in the data set with 11628 observations are verified by experiments. The algorithm has an accuracy of nearly 94% and a 0.96 Matthews correlation coefficient (MCC) score.

Plasmid DNA is common in bacterial populations and is used by synthetic biologists to alter the genetic makeup, and therefore function, of cells. At each cell division in a plasmid containing population, there is a probability of plasmid loss, giving rise to a differentiated population. The plasmid loss rate is difficult to measure, because it is small and quickly overshadowed by exponential growth of the subsequent plasmid free population. In addition, biologists observe complex intracellular plasmid dynamics, involving 1) formation of plasmid clusters that reduce plasmid diffusion, 2) plasmid localization to cell poles, and 3) plasmid replication with negative feedback. Mathematical models are useful in understanding how the plasmid loss rate is determined by these dynamics, but no previous works have incorporated this level of mechanistic detail. We will discuss simulation studies on the influence of the three mechanisms on the probability of plasmid loss, and under what conditions these effects can be captured by tractable mathematical models.

Many types of antibiotics are believed to cause DNA damage in bacteria. The bacterial SOS response is known to promote bacterial survival during antibiotic treatment by inducing the expression of proteins that repair DNA damage. However, the mechanisms by which antibiotics generate DNA damage and trigger the SOS response remain unclear. Here, we propose a delay differential equation model that predicts the temporal dynamics of the SOS response, under action of ciprofloxacin, a DNA-damaging antibiotic. We calibrate the model using ABC-SMC, with data from time-resolved single-molecule and single-cell microscopy experiments. The model allows us to grants insight into how antibiotic treatments induce complex cell behaviour, such as temporal variation and cell-to-cell heterogeneity in the SOS response.

Organoids, mini engineered tissues, have become increasingly popular in recent years1. Indeed, we are experiencing an explosion of interest in organoids as three-dimensional test beds for biological experiments due to their complex structure and ability to mimic in-vivo tissues2. As a result, work must be done to accurately grow and develop organoids. However, the development of organoids, like all biological tissues, is sensitive to the mechanical signals that can influence behaviour from cell growth to determining cell type and shape3. These mechanical cues can even override biochemical signalling in directing type specification of stem cells4. This is made more complex in multicellular structures where mechanical signals operate over multiple length scales. Thus, mathematical models can provide an elegant framework to shed light on the underlying mechanics. One class of models are those founded in a consideration of continuum elasticity5 as applied to soft tissue mechanics, providing the opportunity to investigate the role of key mechanical factors6. The challenge is to produce models that can capture the active behaviour of cells and their ability to generate force as well as to describe the passive mechanical interactions of the system. We present here a model that captures key force generating mechanisms of organoids, namely cell contractility and cell growth. We describe the interaction between contractility and tissue growth and how their antagonistic behaviour can introduce key mechanical signals that may influence behaviour. As a final step, we consider the potential mechanisms by which mechanical feedback into cell control can be incorporated into our model and the impact this will have. References 1. Schutgens, F. & Clevers, H. Human Organoids: Tools for Understanding Biology and Treating Diseases. https://doi.org/10.1146/annurev-pathmechdis-012419-032611 15, 211–234 (2020). 2. Kim, J., Koo, B. K. & Knoblich, J. A. Human organoids: model systems for human biology and medicine. Nature Reviews Molecular Cell Biology vol. 21 571–584 Preprint at https://doi.org/10.1038/s41580-020-0259-3 (2020). 3. Orr, A. W., Helmke, B. P., Blackman, B. R. & Schwartz, M. A. Mechanisms of Mechanotransduction. Dev Cell 10, 11–20 (2006). 4. Engler, A. J., Sen, S., Sweeney, H. L. & Discher, D. E. Matrix Elasticity Directs Stem Cell Lineage Specification. Cell 126, 677–689 (2006). 5. Taber, L. Alan. Nonlinear theory of elasticity applications in biomechanics. Nonlinear theory of elasticity applications in biomechanics (World Scientific). 6. Littlejohns, E. & Dunlop, C. M. Mechanotransduction mechanisms in growing spherically structured tissues. New J Phys 20, 043041 (2018).

This talk addresses the mathematical modelling of phototrophic granular biofilms, spherical, dense aggregates constituted by a relevant phototrophic component and developed in presence of light. These biofilm granules are typically cultivated within bioreactors, and represent an innovative technology in the field of wastewater treatment. Specifically, the presented model describes both the growth of phototrophic granules and the related wastewater treatment process occurring in the bioreactor. The biofilm granule has been modelled as a free boundary domain with radial symmetry, which evolves over time as a result of microbial growth, attachment and detachment processes. Hyperbolic and parabolic partial differential equations (PDEs) have been considered to model at mesoscale the transport and growth of sessile biomass and the diffusion and conversion of soluble substrates. The macroscale behaviour of the system has been modelled through first order impulsive ordinary differential equations (IDEs), which reproduce a sequencing batch reactor (SBR) configuration. Phototrophic biomass has been considered for the first time in granular biofilms, and cyanobacteria and microalgae have been accounted separately, to model their different growth and granulation abilities. To describe the key role of cyanobacteria in the photogranulation process, the attachment velocity of all suspended microbial species has been modelled as a function of the cyanobacteria concentration in suspended form. The model takes into account the main biological processes involved in photogranules-based systems: metabolic activity of cyanobacteria, microalgae, heterotrophic and nitrifying bacteria, microbial decay, EPS secrection, symbiotic and competitive interactions between different species, light-dark cycle, light attenuation across the granule and photoinhibion phenomena. The model has been integrated numerically, and the results show its consistency in describing the photogranules evolution and ecology, and highlight the advantages of the photogranules-based technology, analyzing the effects of the influent wastewater composition and light conditions on the process.

Antimicrobial resistance is a global health concern that requires all possible means of control it. To avoid the onset of multidrug-resistant strains each drug is subject to a maximum time of administration and to a minimum time of administration for effectiveness, referred to as Waiting Time Constraints (WTCs) on biomedical treatments. Treatment with drug combinations and appropriated WT specifications can be modeled by a nonlinear switched system, where a mode (or subsystem) represents the specific administrated drug and the schedule of drugs is associated with an optimal control problem that aims to reduce therapeutic escape. From a dynamic perspective, WTCs significantly alter the regions of the state space of the system that can be feasibly stabilized as they prevent excessive or insufficient time spent in a given mode. Indeed, the literature lacks results on this subject except for recent outcomes for the linear case. Understanding the regions that can be feasibly stabilized is as important as control developing or system modeling; without them no well-posed control can be formulated. Nonetheless, the regions that can be feasibly stabilized by predictive controllers for a nonlinear switched system of antimicrobial resistance, where drugs are subject to WTCs, remain unknown. In this work we present and analyze a series of algorithms to compute stabilizing regions for a switched mathematical model under WTC. The results applied to antibiotic-sensitive and resistance bacteria dynamic population during the course of multi-antibiotic treatment of an infected host addresses the following problems: (i) the existence and characterization of general regions of the state space wherein controlled states trajectories under WTC can feasibly (and indefinitely) remain inside; (ii) when the conditions on WTC to avoid the emergence of resistance allows the existence of feasible and stable control strategies for the success of multiple drug treatment in suppressing the infection; and (iii) when the condition on WTC for the treatment regimen predicts the failure of the treatment due to resistance.

Bacteria can form dense communities called biofilms, where cells are embedded in a self-produced extracellular matrix. Exploiting competitive interactions between strains within the biofilm context can have potential applications in biological, medical, and industrial systems. By combining mathematical modelling with experimental assays, we reveal that spatial structure and competitive dynamics within biofilms are significantly affected by the location and density of the founder cells used to inoculate the biofilm. Using a species-independent theoretical framework describing colony biofilm formation, we show that the observed spatial structure and relative strain biomass in a mature biofilm comprising two isogenic strains can be mapped directly to the geographical distributions of founder cells. Moreover, we define a predictor of competitive outcome that accurately forecasts relative abundance of strains based solely on the founder cells potential for radial expansion - a result we confirmed experimentally. Consequently, we reveal that variability of competitive outcome in biofilms inoculated at low founder density is a natural consequence of the random positioning of founding cells in the inoculum. Extension of our study to non-isogenic strains that interact through local antagonisms, shows that even for strains with different competition strengths, a race for space remains the dominant mode of competition in low founder density biofilms. Our results, verified by experimental assays using Bacillus subtilis, highlight the importance of spatial dynamics on competitive interactions within biofilms and hence to related applications

Understanding the evolution of antimicrobial resistance (AMR) is central for their treatment. In this talk, I want to show a possible way to address this problem from a statistical point of view, namely the hypercubic inference, which we developed and introduced during the last years at the University of Bergen. The basis of this model is a hypercubic transition graph, whose nodes represent possible resistance states and the edges between correspond to the different evolutionary steps. This new approach allows us to efficiently make predictions about the most likely evolutionary pathways leading to AMR and learn their structure and variability, even if we have incomplete datasets with uncertain states. For this we can either use Bayesian inference via Monte Carlo Markov Chain methods or a frequentist approach for the estimation of likelihoods, whereby we only need cross-sectional datasets. While we focus here on AMR, hypercubic inference can be and has been used in a very wide range of problems involving evolutionary accumulation and disease progression, including ovarian cancer, severe malaria, genome and behavioral evolution, and educational progress. The focus of the talk will be the introduction and explanation of the methods themselves, whereby I will address both the advantages and strengths of using a hypercubic structure, but also open problems and ongoing work. In addition, I will also present the results of concrete current applications to real AMR datasets from Klebsiella pneumoniae and Escherichia coli and discuss some biological insights that can be derived from them.

Cyanobacteria are the largest group of photosynthetic organisms on earth. They can survive in very severe conditions (e.g., in deep oceans and near poles) due to complex mechanisms that help them adapt to the specific light spectrum of their surroundings. They are evolved to adjust their metabolism and photosynthesis optimally to their environment. One such genetic switch found in cyanobacteria is the blue-green light switch: a simple system that is likely transferrable to other bacterial species. Thus, this switch has potential use in new regulatory systems for biotechnology and optogenetics. Additionally, cyanobacteria are phototrophs and are already being used to develop more sustainable biotechnology platforms. Characterizing and optimizing a highly responsive gene regulatory system that works efficiently in individuals and populations of cyanobacteria will help to advance their usefulness in biotechnology and in the production of biofuels. We developed an agent-based model of cyanobacteria mixed populations using PhysiCell, an open source physics-based modeling software. We explored the potential of combining opto-genetic and diffusible chemical controls to guide novel spatiotemporal pattern formation. We experimented with different blue-green light spectrums along with population density and cell motility parameters to probe this system for light mediated spatial pattern formation.

The loss of biodiversity due to the spread of destructive tree diseases and invasive pests within forests across the world is having an enormous environmental, economic, and social impact. Enhancing biosecurity is a key priority, through the control of existing diseases and pests, and by building forest resilience against new ones. Thus, we are working in collaboration with multiple forestry partners to develop mathematical models to deepen our understanding of the fundamental behaviours of key pests and pathogens, act as predictive tools for forecasting, and to explore different control strategies. Broadly, we use a combination of partial differential equations, agent-based modelling, and statistical inference techniques. This talk will present an overview of the collaborative work to date, including a case study example of the oak processionary moth infestation in London.

Brassicaceae plants have the glucosinolate-myrosinase defense system, jointly active against herbivory. Glucosinolates (GLS) are hydrolyzed by myrosinase to produce isothiocyanates as soon as herbivory begins. Isothiocyanates exert detrimental effects on the feeding insects. However, constitutive GLS defense is observed to occur at levels that do not deter all insects from feeding. That prompts the question of why Brassicaceae plants have not evolved a high constitutive defense. The answer may lie in the contrasting relationship between plant defense and host plant preference of specialist and generalist herbivores. One of the reasons plants are in this dilemma is that they do not know what kind of herbivore will attack them in any given year, and thus have to be prepared for different possibilities. GLS content increases the susceptibility to specialist insects because these are attracted to plants with a high GLS content and are capable of coping with the toxin. In contrast, generalists are deterred by the plant GLS. Although GLS can attract the natural enemies (predators and parasitoids) of these herbivores, enemies can reduce herbivore pressure to some extent only. So, plants can be overrun by specialists if GLS content is too high, whereas generalists can invade the plants if it is too low. Therefore, an optimal constitutive plant defense can minimize the overall herbivore pressure. To explain optimal defense theoretically, we represent the contrasting host selection behavior of insect herbivores and, in addition, the emergence of their natural enemies by a non-autonomous ordinary differential equation model, where the independent variable is the plant GLS concentration. From the model, we quantify the optimal amount of GLS, which minimizes the total herbivore (specialists and generalists) pressure. That quite successfully explains the evolution of constitutive defense in plants from the perspective of optimality theory.

The spread of American Bullfrog, one of the 100 of the World’s Worst Invasive Alien Species, has a great impact on the surrounding ecosystem, so it will be very important to find out the possibility of spread by region. We assess whether bullfrogs will continue to spread, stop spreading and maintain populations, or become extinct 60 years after their introduction to Korea. This study is based on the results of national surveys that observed the distribution. The entire data is divided into 25 regional clusters using the Hierarchical clustering method, and the degree of spread is predicted by CNN(Convolution Neural Network) method which trains and learns the rules of ECA(elementary cellular automata) that determine evolution of the clusters. We predict the probabilities of the ECA rules for each cluster. The mean value of the population according to the predicted rules is defined as the spreading intensity and evaluated, which is multiplied by the habitat suitability to get an assessment of bullfrog spreading. Habitat suitability is obtained using Maxent.

Hepatitis D virus (HDV) and hepatitis B virus (HBV) coinfection is the most severe form of chronic viral hepatitis. HDV is considered a satellite virus because it relies on hepatitis B surface antigen (HBsAg) to propagate. Treatment against HDV chronic infection with pegylated interferon-α2a (pegIFN) therapy is suboptimal and affects both HDV and HBV. The investigational drug called lonafarnib (LNF) targets HDV only, providing a unique opportunity to study the interplay between HDV and HBV. We recently developed a mathematical model to study the interplay between HDV and HBV in chronic HDV/HBV patients [1]. I will present our efforts to characterize the frequent kinetic data of HDV, HBV, HBsAg, and LNF pharmacokinetics obtained from 15 patients who were treated with LNF, LNF+ritonavir, or LNF+pegIFN [2].In addition, I will present our modeling efforts in extending our model [1] to account also for HBsAg kinetics and to estimate HDV/HBV kinetic parameters and LNF±pegIFN efficacies using both individual and population fitting approaches.

Eastern Equine Encephalitis Virus (EEEV) is an arthropod-borne, single-stranded positive-sense RNA virus that poses a significant threat to public health and national security. Compared to similar viruses such as SARS-CoV-2 or Hepatitis C virus, EEEV causes severe encephalitis when neuroinvasive leading to a human mortality rate of ~30-70%. Moreover, there are no preventative or standardized therapies, leaving patients to rely solely on supportive care. In addition, studies have shown that EEEV is easily aerosolized making it an ideal biowarfare agent. Although the molecular components and interactions of infection, replication, and amplification of EEEV within the host cell are well-studied, how these mechanisms integrate to determine the dynamics of RNA viral replication and host immune responses remains unclear, limiting our ability to advance therapeutic development. Computational models provide a powerful tool for probing both quantitative and qualitative effects arising from the modulation of viral infections. Here, we present a systems modeling approach to elucidate the mechanisms regulating the precise dynamics of EEEV replication through the development of a mechanistic mathematical model. Specifically, this model describes attachment, entry, uncoating, replication, assembly, and export of both infectious virions and virus-like particles within mammalian cells. The model recapitulates known characteristics of EEEV infection, including the timing and amplitude of virion production, and identifies genome replication as the significant rate-limiting step during infection. Additionally, this model highlights the possibility, which will be tested experimentally, that a mismatch between the production of viral RNA and viral proteins could result in the inability to produce infectious virions 12 hours post-infection. We are currently working to expand the model to incorporate type I interferon induction within an infected host cell. This will provide a comprehensive perspective on the conditions required for maximizing host response efficacy and determine the key steps of immune system activation required for successful suppression of viral infection.

After infecting a host, a viral strain may increase rapidly within the body and produce mutants with a faster proliferation rate than the virus itself. However, most of the mutants become extinct because of the stochasticity caused by the small number of infected cells. In addition, the mean growth rate of a mutant strain decreases with time because the number of susceptible target cells is reduced by the wild-type strain. In this study, we calculated the fraction of mutants that can escape stochastic extinction, based on a continuous-time branching process with a time-dependent growth rate. We analyzed two cases differing in the mode of viral transmission: (1) an infected cell transmits the virus through cell-to-cell contact with a susceptible target cell; (2) an infected cell releases numerous free viral particles that subsequently infect susceptible target cells. The chance for a mutant strain to be established decreases with time after infection of the wild-type strain, and it may oscillate before convergence at the stationary value. We then calculated the probability of escaping stochastic extinction for a drug-resistant mutant when a patient received an antiviral drug that suppressed the original strain. Combining the rate of mutant production from the original strain and the chance of escaping stochastic extinction, the number of emerging drug-resistant mutations may have two peaks: one soon after the infection of the original type and the second at the start of antiviral drug administration. Hayashi R, Iwami S, and *Iwasa Y. 2022. Escaping stochastic extinction of mutant virus: temporal pattern of emergence of drug resistance within a host. Journal of Theoretical Biology 537, 111029. Hayashi R., and Iwasa, Y. Temporal pattern of the emergence of a mutant virus escaping cross-immunity and stochastic extinction within a host. (in review)

Coral reefs, among the most diverse ecosystems in the ocean, currently face major threats due to multiple stressors such as pollution, unsustainable fishing practices, and perturbations in environmental parameters brought on by climate change. Reefs also sustain regular wounding from other sea life and human activity. Recent reef preservation practices have even involved intentional wounding by systematically breaking coral fragments and relocating them to revitalize damaged reefs. Despite its importance, very little research has explored the inner mechanisms of wound healing in corals. Some reef-building corals have been observed to initiate an immunological response similar to those observed in humans and other vertebrates. Utilizing past models of inflammation and early proliferation and remodeling, we formulated a mechanistic model for wound healing in corals. The model consists of four differential equations mediating wound debris, inflammation, proliferation, and wound closure. The model is coupled with experimental data for linear and circular shaped wounds on Pocillopora damicornis fragments. A preliminary parameter set was obtained by fitting to the wound closure times obtained empirically and to expected temporal trends observed in other coral species and in humans and other vertebrates. A variety of mathematical methods were applied for model analysis including local sensitivity analysis. Results were used to define an identifiable set of parameters. The parameter space was also explored to exhibit the diverse model outcomes and their biological implications. Keywords: stony corals, inflammation, differential equations, parameter estimation, sensitivity analysis

The pandemic preparedness of nursing homes has been a major concern for decades. The COVID-19 pandemic proved that the concerns were valid, as it caused devasting death tolls in nursing home facilities. This study presents an agent-based modeling framework to better understand the dynamics of pandemics within and between nursing homes. This is sharply distinct from many agent-based modeling works that resemble the spread of the infection within a single nursing home. We first calibrate the model of multiple nursing homes using the available COVID-19 data. Then we investigate the effects of shared staff on the efficacy of Covid-19 preventive policies through extensive simulations. It is shown that shared staffing can significantly diminish the efficacy of preventive policies. In conclusion, the nursing workforce is a determining factor for pandemic preparedness.

Background: The COVID-19 pandemic has become a new global public health crisis, and to large extent, its capacity to cross natural geographic barriers is attributed to human mobility and contact patterns which vary with time and specific locations. Therefore, an agent-based model (ABM) which relates populations mobility patterns in different locations in compliance with on site COVID-19 control measures is proposed to investigate how opening and closing protocols would have been best implemented in Mozambique. Methods: For spatial dynamics, a survey was carried out in the city of Maputo as a case study to estimate populations mobility patterns and contact matrices among individuals in different locations (home, school, work place, worship place, market and any other place of gathering) during specific periods of the day (morning, afternoon and night) for both week days and weekends. Individuals are explicitly represented by agents associated to disease characteristics and their decision to remain or move to a new place is based on a probability estimated from the survey and on site declared control measures. Results: The results show that at $50%$ of social distancing compliance, complete lockdown of schools, workplaces, worship places with exception of markets is the only scenario that result in the reduction and shift of the peak by $3%$ and 3 days respectively. School closure showed significant effect that at $75%$ and $85%$ of social distancing adherence resulted in the reduction and shift of the peak by $15%$ and 4 days, and $51%$ and 24 days respectively. While closure of worship places rendered little effect due to limited frequency and duration of activities in a given location. Conclusions: This study has demonstrated the use of simulation models to investigate the implementation of opening and closing policies for the control of COVID-19 pandemic at local scale by leveraging between the mobility of individuals and adherence to social distancing.

Percent positivity, the ratio of positive SARS-CoV-2 tests to total number of tests, has been used throughout the COVID-19 pandemic as a proxy for the current level of transmission in a community. Simultaneously, wastewater SARS-CoV-2 monitoring has been implemented, but is a highly variable metric whose direct utility has yet to be fully explored. As we transition from pandemic response to endemic management, testing efforts have been reduced and the predictive value of test percent positivity has been called into question. We build a series of models incorporating SARS-CoV-2 test positivity, wastewater SARS-CoV-2 levels, and syndromic surveillance data streams to explore changing transmission dynamics. A county level model is developed to forecast hospitalizations and tested against an ARIMA based on hospitalizations alone. A 21-day forecast is developed with sliding scale cross validation. We validate and quantify uncertainty in commonly used public health metrics and explore differences in model selection between variants. Data from the winter 2022-23 season are reserved as a final test for the model. In this work, we examine how to effectively predict hospitalizations in a changing testing environment.

We investigate an SIRS epidemic model with spatial diffusion and nonlinear incidence rates. We show that for small diffusion rate of the infected class DI , the infected population tends to be highly localized at certain points inside the domain, forming K spikes. We then study three distinct destabilization mechanisms, as well as a transition from localized spikes to plateau solutions. Two of the instabilities are due to coarsening (spike death) and self-replication (spike birth), and have well-known analogues in other reaction-diffusion systems such as the Schnakenberg model. The third transition is when a single spike becomes unstable and moves to the boundary. This happens when the diffusion of the recovered class, DR becomes sufficiently small. In all cases, the stability thresholds are computed asymptotically and are verified by numerical experiments. We also show that the spike solution can transit into an plateau-type solution when the diffusion rates of recovered and susceptible class are sufficiently small. Implications for disease spread and control through quarantine are discussed.

We study a diffusive epidemic model and examine the spatial spreading dynamics of a multi-strain infectious disease. In particular, we address the questions of competition-exclusion or coexistence of the disease's strains. Our results indicate that if one strain has its local reproduction function spatially homogeneous, which either strictly minimizes or maximizes the basic reproduction numbers, then the phenomenon of competition-exclusion occurs. However, if all the local reproduction functions are spatially heterogeneous, all strains may coexist. In this case, we provide complete information on the large time behavior of classical solutions for the two-strain model when the diffusion rate is uniform within the population and the ratio of the local transmission rates is constant. Particularly, we prove the existence of two critical superimposed functions that serve as threshold values for the ratio of the transmission rates and that of the recovery rates. Furthermore, when the populations' diffusion rates are small, our results on the asymptotic profiles of coexistence endemic equilibria indicate a spatial segregation of infected populations.

One-dimensional discrete-time population models, such as Logistic or Ricker growth, can exhibit periodic and chaotic dynamics. Incorporating epidemiological interactions through the addition of an infectious class causes an interesting complexity of new behaviors. Previous work showed that infection that abrogates fecundity can lead to unexpected increases in total population size, a phenomenon known as the ‘hydra effect.’ Here, we examine a two-dimensional susceptible-infectious (SI) model with underlying Ricker population growth and show that the disease-free system has a distinct bifurcation structure from the system with infection. We use numerical bifurcation analysis to determine the influence of infection on the types and appearance of qualitatively distinct long-time dynamics. We find that disease-induced mortality leads to the appearance of multistability, such as stable four-cycles and chaos dependent upon the initial condition. In addition, we examine the appearance and extent of the hydra effect, particularly when infection is introduced during cyclic or chaotic population dynamics.

Aedes aegypti and Aedes albopictus mosquitoes are the causative agents of dengue (DENV), chikungunya (CHIKV), and Zika (ZIKV) virus infections in humans. The co-circulation of at least two viruses/serotypes, which is common in countries worldwide, such as Columbia and Brazil in Latin America, can cause potential interactions among the viruses/serotypes and misdiagnosis in the lack of adequate laboratory tests due to similar clinical symptoms among their disease courses. We generalized a deterministic compartmental model to analyze how each disease dynamics changes under the potential antagonistic or synergistic interaction among the viruses/serotypes. Our simulation studies showed that under no DENV vaccine, vector control, and interaction among the viruses/serotypes, the peak of the incidence rates for people with no prior infections happens earlier than those cases with one or two prior infections, mostly because the proportion of fully susceptible people is larger than people with at least one prior infection. We observed higher incidence rates for single/multi infections and an earlier peak of the epidemics for single infections when a prior infection by a virus such as ZIKA causes synergistic cross-immunity against CHIKV and DENV serotypes, compared to the situation when it causes antagonistic cross-immunity. Identification of the cross-immunity is not possible when susceptibility statuses of the population are unknown because the high/low incidence rates could be either the results of high/low baseline transmission rates or antagonistic/synergistic interaction effects among the viruses.

Matrix decomposition methods such Principal Components Analysis (PCA) and Non-Negative Matrix Factorization (NMF), as well as non-linear analogues such as t-SNE and UMAP, have become increasingly popular in bioinformatics to perform data dimension reduction. For data that includes a time-course, a more natural representation is as a tensor (a multi-index array). We show that one can reframe the tensor decomposition of sample-by-feature-by-time-course data as a tensor dictionary learning problem, which effectively models each subject as a sum of rank-one matrices we term 'Feature Canonical Trajectories' (FCTs). The benefits of the FCT representation is that it provides not only an embedding and clustering of subjects, but also a model for the data by representing subject gene expression data as a linear combination of canonical trajectories of feature-sets. We show that by reframing the decomposition as a data model we can also identify new metrics to choose a decomposition algorithm and approximation for improved interpretability, and provide an example of this in action by identifying a novel age-specific FCT associated with platelet vaccination response data.

Posttraumatic stress disorder (PTSD) is a debilitating condition triggered by traumatic events. Notable symptom differences exist between combat-exposed veterans, active-duty personnel, and civilian PTSD cases. However, the underlying biological mechanisms remain elusive. This study aims to uncover the shared biological core modules associated with PTSD by leveraging extensive omics data among various PTSD subtypes. To achieve this, we employed the 'COre Module Biomarker Identification with Network ExploRation' (COMBINER) approach on DNA methylation data to identify key network modules of epigenetic modification across PTSD subtypes. These findings not only enhance our knowledge of PTSD's diverse symptomatology but also pave the way for the development of biomarkers and personalized treatments.

In recent years, a growing interest in personalized medicine has emerged as a result of significant advancements in the fields of biology, data science, and computational modeling. One emerging concept that has obtained attention from the scientific community is the patient digital twin (DT), which aims to develop a comprehensive model to enable clinicians to systematically analyze the complexity of each patient, simulate treatment outcomes, and select the optimal treatment option. Constructing a patient digital twin (DT) involves creating a detailed computational model that can capture the unique characteristics of an individual patient. The model should include information about the patient's medical history, current health status, genetic makeup, and other relevant factors that can affect treatment outcomes. However, obtaining all this information can be challenging, as clinical data on individual patients is often limited. To overcome this limitation, researchers can use mechanistic models to replicate observed phenomena across various scenarios. In this study, we developed a mechanistic model of cancer-immune interactions in pulmonary micrometastasis. We found that the model could express a wide range of patient trajectories, from complete tumor elimination to uncontrollable growth. Using high-throughput computing platforms, we analyzed 100k virtual patient trajectories by exploring the parameter space of this model. Initially, we analyzed patient data at a single time point using clustering methods, but the results did not clearly identify patient trajectories. Further investigation revealed that the same patient could have completely different trajectories, making it challenging to categorize patients. The mechanistic model helped us understand this issue, showing that early or non-interactions between macrophages and invading tumor cells were responsible. This also explained the limitations of traditional patient stratification based on data alone. Additionally, it highlighted the need for digital twins that are patient-specific, dynamical, and continuously updated with new patient data instead of one-time calibration.

Single-cell RNA sequencing (scRNA-seq) has enabled various analyses, including cellular phenotyping and gene regulatory network reconstruction. However, the sparsity, high dimensionality, bias, skewness in data distribution, and technical noise in scRNA-seq data present challenges for downstream analyses. In order to tackle these issues, conventional packages depend on log-normalization for preprocessing and require users to select the reduced dimension when employing various dimensionality reduction methods. Nonetheless, these approaches can result in signal distortion and subjectivity when determining reduced data dimensions. To overcome the limitation of conventional methods, we developed a new tool that corrects signal distortion during preprocessing and effectively filters out noise in data, enhancing the reliability of the outcome of dimensionality reduction. Our tool demonstrated superior performance compared to 9 widely used packages, including Seurat, Scanpy, and Monocle3 when tested on 53 real and simulated datasets.

In volume transmission (or neuromodulation) neurons do not make one-to-one connections to other neurons, but instead simply release neurotransmitter into the extracellular space from numerous varicosities. Many well-known neurotransmitters including serotonin (5HT), dopamine (DA), histamine (HA), Gamma-Aminobutyric Acid (GABA) and acetylcholine (ACh) participate in volume transmission. Typically, the cell bodies are in one volume and the axons project to a distant volume in the brain releasing the neurotransmitter there. We introduce volume transmission and describe mathematically two natural homeostatic mechanisms. In some brain regions several neurotransmitters in the extracellular space affect each others' release. We investigate the dynamics created by this comodulation in two different cases: serotonin and histamine; and the comodulation of 4 neurotransmitters in the striatum and we compare to experimental data. This kind of comodulation poses new dynamical questions as well as the question of how these biochemical networks influence the electrophysiological networks in the brain.

Two ubiquitous components, coupling and noise, may drive complex neural networks to exhibit emergence dynamics. While the roles of equal coupling and symmetric noise have been extensively studied, the general mechanisms of unequal coupling strength and asymmetric noise remain unclear. In this work, we investigate the simultaneous interplay of unequal coupling and asymmetric noise in the simplest network motif of two bi-directionally coupled neural oscillators, each with its own intrinsic noise. Our findings show that noise-induced synchrony can be maximized when one oscillator (source) with weak intrinsic noise is strongly connected to the other oscillator (receiver) with strong intrinsic noise. Furthermore, we extend our study to three coupled neural oscillators with a feed-forward-loop schematic. These results shed new light on the complex interplay between coupling and noise in neural networks.

Gonadotropin-releasing hormone (GnRH) neurons work as a trigger of the reproductive axis in mammals. These neurons exhibit two types of exocytosis: a surge and a pulsatile one. Traditionally, changes in the neuron dynamics are connected to and explained by changes in parameters of the action potential, transmitted by a neuron's membrane. However, in case of GnRH neurons, the experimental data demonstrates that the switch in the type of the hormone release is determined rather by the location of the GnRH neuron activation. Action potential initiated in the proximity of soma is necessary for the surge of GnRH. The second type, the pulsatile release of GnRH, is driven by the synaptic activities on the distal part of the neurons. Both types of the exocytosis initiation target the intracellular calcium dynamics. The increase in calcium ions due to the electrical spikes near soma is short-lived. On the other hand, the increase in calcium ions in the distal parts of the GnRH neurons lasts for tens of minutes. We have built the mathematical and computational models of the electrical and chemical dynamics in GnRH neurons. The model, in silico, reveals the connection between the action potential, neuropeptides, and calcium ion dynamics. In addition, our model confirms the functionality of the bundling between multiple GnRH neurons and its effect on exocytosis synchronization.

Modeling the human cortex is challenging due to its structural and dynamic complexity. Spiking neuron models can incorporate many details of cortical circuits but are computationally costly and difficult to scale up, limiting their scope to small patches of cortex and restricting the range of phenomena that can be studied. Alternatively, one can use simpler phenomenological models, which are easier to build and run but are more difficult to compare directly to experimental data. This talk presents an efficient modeling strategy that aims to strike a balance between biological realism and computational efficiency. The proposed modeling strategy combines a coarse-grained representation with local circuit dynamics to compute the steady-state cortical response to external stimuli. A crucial observation is that as a consequence of anatomical structures and the nature of neuronal interactions, potential local responses can be computed independently of dynamics on the coarse-grained level. We first precompute a library of steady-state local responses driven by possible lateral and external input. Then, the fixed point of the coarse-grained model can be computed by an iterative scheme combined with fast library lookup. Our approach is tested on a model of primate primary visual cortex (V1) and successfully captures essential V1 features such as orientation selectivity. Time permitting, I will also discuss a related project in which we devised an efficient way to explore the parameter space of a primate V1 model, identifying the set of viable parameters as a 'thickened' codimension-1 submanifold of parameter space.

In recent years, game theory has become a valuable tool for studying the evolutionary dynamics of cancer and their response to treatment. In this study, we use an evolutionary game theory model to investigate the interplay between treatment and the environment in the development of drug resistance in tumor populations. The model consists of two tumor populations: a sensitive and a resistant phenotype, and a stromal population that interacts with the tumor cells. We use replicator equations to simulate the evolutionary dynamics of the populations and analyze the effects of different treatment strategies on the development of drug resistance. To validate the model, we integrate experimental data from previous studies that have investigated the role of stromal cells, specifically carcinoma-associated fibroblasts, in weakening the tumor's response to treatment in cancer. We use simulations and analyses to optimize adaptive therapies in this context by adjusting treatment frequency and dosage. Our results demonstrate that adaptive therapies using a combination of drugs can significantly reduce the emergence of drug resistance in the tumor population. Furthermore, our simulations suggest that adjusting the treatment frequency and dosage can further optimize the adaptive therapies, taking into account the interactions between the tumor and stromal cells. These results highlight the importance of considering the evolutionary dynamics of cancer cells when designing treatments. By understanding how selective pressures affect tumor populations and how they respond to treatment, we can develop more effective therapies that can prolong patient survival and improve their quality of life. Moreover, the evolutionary game theory model presented here can be adapted to other cancer types, potentially leading to the development of more personalized and effective treatments.

Bone metastasis is a lethal consequence for prostate cancer patients, mostly due to the emergence of resistance and therapy failure. Bone-targeting radiotherapy with Radium223 (Rad223), a radioisotope emitting genotoxic alpha-radiation with limited tissue penetrance (∼100 µm), prolongs the survival of patients with metastatic prostate cancer (PCa). The clinical response to Rad223 is often followed by detrimental relapse and progression. Whether Rad223 causes tumor-cell directed cytotoxicity in vivo remains unclear, and, additionally, effective strategies to improve long-term Rad223 efficacy have not been developed yet. Integrins are heterodimeric transmembrane receptors that, through their activity, support cell growth, decrease cell death and enable radio-resistance mechanisms on exposure to ionizing radiation. Anti-β1 integrin (β1I) targeting improves irradiation treatment outcomes in breast cancer cells and other subcutaneous xenografts. We hypothesize that: i) limited radiation penetrance in situ defines outcome, and ii) targeting β1I would improve 223Ra outcome. In vivo approaches have been limited by huge resource demanding, ethical concerns, and paucity of investigation time points during follow-up. The integration of central biological findings with mathematical modeling allows generating in silico pathophysiological profiles suitable for testing preclinically relevant hypotheses, including predictions on the impact of combinatorial treatments on disease progression. Such predictions can guide preclinical and clinical studies towards more successful outcomes and maximize efficacy. Accordingly, we developed an agent-based model of prostate cancer bone metastasis establishment, growth, and response to putative therapies to: i) predict Rad223 effectiveness in lesions of different sizes, ii) identify Rad223 resistance niches, and iii) optimize Rad223- β1I combinatorial regimen. Our model is regulated via Montecarlo stochastic simulations and seeded with bone metastasis cells endowed with mitosis/apoptosis probability densities which drive the model towards different outcome according to the therapeutic regimen tested. Our in silico model predictions were validated in vivo on humanized bone metastasis mouse model Rad223 delayed the growth of tumors (PC3 and C4-2B cell line) in bone. Cancer cell lethality in response to Rad223 was profound but zonally confined along the bone interface compared with the more distant tumor core, which remained unperturbed. In silico simulations predicted greater efficacy of Rad223 on single-cell lesions and minimal effects on larger, as further confirmed in vivo for PC3 and C4-2B tumors. Micro-tumors showed severe growth delay or eradication in response to Rad223, whereas macro-tumors persisted and expanded. The relative inefficacy in controlling large tumors points to application of Rad223 in secondary prevention of early bone-metastatic disease and regimens co-targeting the tumor core. Interference with β1I combined with Rad223 reduced PC3 cell growth in bone and significantly improved overall mouse survival, whereas no change was achieved in C4-2B tumors. Anti-β1I treatment decreased the PC3 tumor cell mitosis index and spatially expanded Rad223 lethal effects 2-fold, in vivo and in silico. Regression was paralleled by decreased expression of radio-resistance mediator. Targeting β1I significantly improves Rad223 outcome and points toward combinatorial application in PCa tumors with high β1I expression.

Tumor relapse during therapies is thought to reflect the ability of tumor cells to escape treatment via cell-intrinsic genetic and epigenetic changes, either preexistent or evolved de novo under the drug-imposed cytostatic and cytotoxic selective pressures. Multiple published and unpublished reports, show that at least some of the otherwise drug-sensitive tumor cells might be able to avoid elimination due to microenvironmental factors, such as pro-survival paracrine signals, capable of providing alternative means of survival despite the shutdown of oncogenic signaling. Concurrently, clinical studies indicated that higher stroma-to-tumor content has a poorer prognosis and increased risk of relapse in various cancers. However, despite our advances in deciphering the molecular mechanisms behind tumor progression, our understanding of how stroma contributes to the evolution of tumor resistance and relapse remains limited. Our objective is to build stochastic mathematical models to help understand the stroma-tumor eco-evolutionary interplay so that we can design new therapeutic strategies for tumor elimination or treatments that provide long-term control of tumor growth. The models integrate the spatiotemporal stroma-tumor cells interactions measured using a novel purposely-designed spatial analysis pipeline and combine it with in vivo and in vitro tumor growth dynamics data. We focus on two distinct stromal effects studied in our lab in vitro and in vivo, namely drug-sheltering effects against targeted therapies in Non-Small Cell Lung Cancer (NSCLC) and enhanced tumor cells’ proliferation, indiscriminately of treatment presence, in Triple Negative Breast Cancer (TNBC). The preliminary results show that enhanced proliferation in the vicinity of stroma in TNBC could be sufficient to drive tumor relapse after four cycles of chemotherapy and greatly accelerates relapse in tumors where the treatment has a low killing rate. In the context of NSCLC, the drug-sheltering effects of stroma could lead to tumor relapse in absence of hard-wired resistance, and disrupting the drug-sheltering stromal effects result in a higher rate of tumor eradication. Our results show that by using carefully designed and calibrated mathematical models, we can gain a deeper understanding of the ecological mechanisms that lead to tumor relapse as well as uncover new therapeutic strategies that account for stromal effects and are successful at eradicating tumors.

Engineered T cell receptor (TCR) T cells are expected to drive strong anti-tumor responses upon recognition of the specific tumor antigen, with rapid expansion and cytotoxic functions, causing tumor cell death. However, although TCR T cell therapy against cancers is promising, it remains difficult to predict which patient will have better therapeutic outcome and why. We develop a mathematical model to re-assess some mechanisms of insufficient efficacy of TCR T cell therapy in HPV+ cervical cancer model of mice. We consider a dynamical system that follows the population of cancer cells, effector TCR T cells, regulatory T cells (Tregs), and 'non-tumor killing' TCR T cells. We demonstrate: 1) the majority of TCR T cell within the tumor is 'non-tumor killing' TCR T cells, such as exhausted cells, which are highly active but contribute little or have no direct cytotoxicity in the tumor microenvironment (TME); 2) there are two important conditions for tumor regression: the reversal of the immunosuppressive TME by depleting Tregs, and the increased proliferation of effector TCR T cells with antitumor activity. Using mathematical modeling, we show that certain treatment protocols have the potential to improve therapy responses.