ABSTRACT. Multiple-objective optimization is common in biological systems. Some of these objectives are incompatible, thus tradeoffs are necessary. living organisms are constantly under selection pressure to maximize their fitness to the environment through optimizing multiple objectives such as growth rate and resistance to environmental fluctuations. In this talk I will focus on two problems that we worked on recently.

In the mammalian olfactory system, each sensory neuron stochastically expresses one and only one out of up to thousands of olfactory receptor (OR) gene alleles; at the organism level, the types of expressed ORs need to be maximized. This Nobel-prize winning observation raises one of the most intriguing puzzles in neurobiology that remains elusive after several decades of intensive investigations: how can both monoallelic and diverse expression of OR be ensured at the same time? Our theoretical studies unraveled how cells achieve these objectives through simple physical principles. The model makes extensive testable predictions. Some are so counterintuitive that the corresponding experimental results, published before the model, received great skepticism in the field, and are predicted by the model with mechanistic explanation.

Acute Kidney Injury (AKI) affects ~13.3 million patients, and is associated with ~1.7 million death globally each year. Even recovered from AKI patients have much higher risk to develop chronic kidney diseases (CKD). About 13% of the population in US suffer from CKD, and the disease eventually leads to end-stage renal disease (ESRD). There is no effective treatment for either AKI or CKD. AKI and CKD are traditionally two separate fields of study. Using combined mathematical modeling and mouse model studies, we showed that the transition from AKI to CKD comes from evolutionary compromise, and resolved some decades-long debate. Furthermore, our studies showed that targeting the recovering dynamics after AKI can potentially improve the repair process without increasing the risk of CKD.

ABSTRACT. While many tumors initially respond to chemotherapy, regrowth of surviving cells compromises treatment efficacy in the long-term. The cell biological basis of this regrowth is not understood. Here, we characterize the response of individual, patient-derived neuroblastoma cells driven by the prominent oncogene MYC to the first-line chemotherapy, doxorubicin. Combining live-cell imaging, cell-cycle-resolved transcriptomics, and mathematical modeling, we demonstrate that a cell’s treatment response is dictated by its expression level of MYC and its cell-cycle position prior to treatment. All Low-MYC cells enter therapy-induced senescence. High-MYC cells, by contrast, disable their bistable cell-cycle checkpoints, forcing renewed proliferation despite treatment-induced DNA damage. After treatment, the viability of High-MYC cells depends on their cell cycle position during treatment: newborn cells promptly halt in G₁ phase, repair DNA damage, and form regrowing clones; all other cells show protracted DNA repair and ultimately die. These findings demonstrate that fast-proliferating tumor cells may resist cytotoxic treatment non-genetically, by arresting within a favorable window of the cell cycle. 

ABSTRACT. Neural stem cells (NSC’s) often give rise to mixed population of cells during differentiation. However, the dynamical origin of these mixed states is poorly understood. In this article, our mathematical modeling study demonstrates that the bone morphogenetic protein 2 (BMP2) driven differential differentiation dynamics of NSC’s in central and peripheral nervous systems essentially function through two distinct bi-stable switches that are mutually interconnected. Stochastic simulations of the model reveal that the mixed population originates due to the existence of these bi-stable switching regulations while the maintenance of such mixed states depends on the level of stochastic fluctuations of the system. Importantly, the model predicts that by individually altering the expression level of key regulatory proteins, the NSC’s can be converted entirely to a preferred phenotype for BMP2 doses that earlier resulted into mixed population. Our findings show that efficient neuronal regeneration can be achieved by systematically maneuvering the differentiation dynamics.

ABSTRACT. The initial rostrocaudal patterning of the neural tube leads to differential expression of Hox genes that contribute to the specification of motor neuron (MN) subtype identity. Although several Hox mRNAs are expressed in progenitors in a noisy manner, these Hox proteins are not expressed in the progenitors and only become detectable in postmitotic MNs. MicroRNA biogenesis impairment leads to precocious expression and propagates the noise of Hoxa5 at the protein level, resulting in an imprecise Hoxa5-Hoxc8 boundary. Using in silico simulations, we uncovered two feed-forward Hox-miRNA loops accounting for the precocious and noisy Hoxa5 expression, as well as an ill-defined boundary phenotype in Dicer mutants. In addition, we identified mir-27 as a major regulator coordinating the temporal delay and the spatial boundary of Hox protein expression. Using more detailed computational analysis, we further predict that the Hox-miRNA circuit underlies a bistable switch with delayed Hoxa5 expression, and that the bistability and the delay synergistically enables a beneficial effect of Hoxa5 transcriptional noise in sharpening of the Hoxa5-Hoxc8 boundary. This intracellular noise contributes to the attenuation of the fluctuations in extracellular morphogens. Therefore, the microRNA can both control and take advantage of the Hoxa5 expression noise. Our results demonstrate a novel mechanism for Hox-miRNA circuit to confer robustness to both individual MN identities and the tissue boundary.

ABSTRACT. Cell polarization toward an attractant is influenced by both physical and chemical factors. Most existing mathematical models are based on reaction-diffusion systems and only focus on the chemical process occurring during cell polarization. However, membrane tension has been shown to act as a long-range inhibitor of cell polarization. Here, we present a cell polar- ization model incorporating the interplay between Rac GTPase, filamentous actin (F-actin), and cell membrane tension. We further test the predictions of this model by performing sin- gle cell measurements of the spontaneous polarization of cancer stem cells (CSCs) and non-stem cancer cells (NSCCs), as the former have lower cell membrane tension. Based on both our model and the experimental results, cell polarization is more sensitive to stimuli under low membrane tension, and high membrane tension improves the robustness and stability of cell polarization such that polarization persists under random perturbations. Fur- thermore, our simulations are the first to recapitulate the experimental results described by Houk et al., revealing that aspiration (elevation of tension) and release (reduction of tension) result in a decrease in and recovery of the activity of Rac-GTP, respectively, and that the relaxation of tension induces new polarity of the cell body when a cell with the pseudopod- neck-body morphology is severed.

ABSTRACT. Rho-GTPases are master regulators of polarity establishment and cell morphology in many eukaryotes. Upon receiving relevant signals, Rho-GTPases become concentrated in clusters at the cell cortex, from where they regulate the cytoskeleton to influence cell behavior. The biochemical mechanisms underlying such clustering include common features such as positive feedback and differential mobility of membrane-associated and cytoplasmic components. The specific functionalities of different cell types require the generation of either one (e.g. the front of a migrating cell) or several clusters (e.g. multiple dendrites of a neuron), but the mechanistic basis for uni-polar or multi-polar outcomes is unclear. Insights into the design principles of Rho-GTPase circuits were provided by simple two-variable reaction-diffusion models that capture essential features of GTPase biochemistry. Here, we use such models to show that when more than one GTPase cluster forms, the core polarity circuit enforces competition between the clusters to yield a uni-polar outcome. However, the efficiency of competition is determined by parameter values in a switch-like manner, with some parameter choices dampening competition to the point that multiple domains can persist on biologically relevant timescales. We derive a “saturation rule” that governs the timescale of competition, and hence whether the system will generate uni-polar or multi-polar outcomes. Our theory suggests that the saturation rule is an fundamental property of the Rho-GTPase polarity machinery, regardless of the specific feedback mechanism.

ABSTRACT. Biological networks are complicatedly wired and therefore we often need more than one control targets to change their attractor state into a desired one. Although various control target selection methods, such as feedback vertex sets (FVS) or minimum dominating sets (MDS) were suggested, the resulting control targets do not indicate synergistic drug targets. The main reason is that such methods are using only the information of network topology. It is true that the steady state or the controllability of a network heavily depends on the network topology, but the synergistic effects are primarily caused by complex network dynamics which are not solely determined by the network topology. Thus, we need to develop a new control strategy that can identify useful synergistic control target pairs based on both topology and dynamics of the network. For this purpose, we have developed a novel method of identifying synergistic control targets by merging the state transition maps before and after virtual perturbations of the network nodes. We also developed a scoring algorithm to evaluate the synergistic effect of each node perturbation. The proposed method compares the average activities and state alteration numbers of network nodes between the desired direction and the opposite direction of the state transition flow in the merged transition map. The scores are weighted based on the phenotypic information according to the attractor classification criteria of the network. We applied the proposed method to published Boolean network models of biological networks and confirmed that our method can identify synergistic drug targets. In addition, by visualizing the merged state transition map and analyzing the state transition flow upon it, we can further reveal the hidden mechanisms of the identified synergistic drug pairs. Acknowledgements: This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Korea Government, the Ministry of Science, ICT & Future Planning (2015M3A9A7067220, 2014R1A2A1A10052404, and 2013M3A9A7046303). It was also supported by the KAIST Grand Challenge 30 Project grant.

ABSTRACT. Cell size in multicellular organisms can vary across several orders of magnitude between different cell types. Even in the unicellular budding yeast, cell volume ranges between 10-200µm3, depending on developmental stage, environmental conditions and age. However, in a given cell type or a specific environment all cells are of the same size. This indicates that size is critical for cell function. How cell size influences cell physiology is not clear.

We used budding yeast to address this question. Using a reversible cell cycle arrest we generated largely oversized cells. Consistent with previous observations in S. pombe, mRNA and total protein levels do not scale with cell volume once cells exceed a critical size, the cytoplasm of oversized cells is therefore diluted.

Increased cell size interfered with basic cellular processes such as cell surface receptor signaling and transcription induction. In addition oversized cells proliferated poorly when released from the cell cycle block. Increasing ploidy or reducing proteasome activity restored size associated signaling-, transcription- and proliferation defects in oversized cells. This shows that dilution of the cytoplasm is detrimental for cell function. We are currently generating in silico models of cell cycle progression in differently sized cells to understand how dilution of critical components affects timing and fidelity of a complex process such as cell cycle progression.

ABSTRACT. Coupled oscillators generate diverse behaviors in a variety of organisms. In autonomous cell systems such as fibroblasts, the circadian clock and the cell cycle are coupled intracellularly with a 1:1 ratio. However, the coupling of clock and mitosis is putatively more complex in heterogeneous, multicellular systems and tissues. Here, we demonstrate dynamic Wnt-mediated intercellular coupling between cell cycle and circadian clock in primary 3D cultures of murine intestinal organoids (enteroids), which is a complex organotypic structures containing intestinal stem cells (ISC), progenitor cells (PC), and differentiated cells (DC). Remarkably, populations of enteroids show circadian clock-dependent synchronized cell division cycles with a period of ~12-h. In contrast, cell cycle measurements from single cells demonstrate a heterogeneous, multimodal distribution of cell cycle times (CCT) with an average period of ~19-h, apparently inconsistent with the population data. To resolve this discrepancy, we developed mathematical model assumed that enteroids consist of ISCs and PCs as proliferating cells with CCTs of 16 and 26-h, respectively, in the absence of circadian connections. With stochastic simulations, we find that circadian rhythms regulate the timing of cell divisions in a heterogeneous population that collectively emerge as synchronized 12-h cell division cycles. Furthermore, when these cells are segregated into ISCs and PCs, we uncovered 1:1 and 1:2 coupling ratios in ISCs and PCs, respectively, in both simulation and experiment. We further observe lack of circadian oscillations in ISCs and PCs, indicating an intercellular signal from DCs mediates circadian clock-dependent synchronized cell division cycles. Simulation and experimental results indicate a key role of circadian rhythms in regulating synchronized divisions of ISCs and PCs via intercellular Wnt signaling.

ABSTRACT. Adaptive immune responses to infection or cancer rely on coordinated programs of cell proliferation and differentiation. Upon infection, naive, antigen-specific T cells expand vigorously and give rise to short-lived effector and long-lived memory cells. Conflicting models have been proposed that suggest either of these subsets to be a precursor of the other; how this subset diversification is regulated by external stimuli like T cell receptor (TCR) avidity, antigen availability or inflammation is largely unknown. Here we show that single cell fate mapping data, interrogated by stochastic population modeling and large-scale model discrimination, are surprisingly informative on both the topology and regulation of differentiation pathways. We find that the phenotypic diversity of T cell responses is generated through stochastic linear cell-fate progression: Naive T cells give rise to slowly dividing memory precursor cells from which rapidly dividing short-lived subsets emerge. This process is modulated but not determined by TCR avidity, which we find to affect the probability with which stochastic division and differentiation events occur. For polyclonal T cell responses, our mathematical model provides a mechanistic explanation for the longstanding observation that high avidity T cell clones within a population of responding T cells only become dominant throughout repetitive immunizations. Proliferation of the T cells is furthermore strongly dependent on both inflammatory signals and continuous stimulation of the TCR. However, we find that the expansion of (central) memory precursors is more dependent on TCR stimuli than the other subsets. Taken together, our mathematical model begins to provide a quantitative picture of the developmental program of T cells during an immune response. Improvements in the quantitative understanding of this process will have implications for immunotherapy and the design of effective vaccines.

ABSTRACT. Proliferating cells tie cell division to growth in order to maintain their size within an optimal range. In budding yeast (S. cerevisiae), size control occurs at START, the point of irreversible commitment to the cell cycle, which can only be passed once a certain critical size is reached. However, how exactly cells measure their size and relay this information to the cell cycle remains controversial.

Here, we present two mechanistically based, mathematical models of size control in yeast cells, in order to assess a recent suggestion that cell growth controls entry into the cell cycle by diluting an inhibitor of START, Whi5. We show that this ‘inhibitor dilution’ model is consistent with most experimental observations but critically fails to account for the size of diploid cells that harbour only one copy of the WHI5 gene. We then propose an alternative model, where Whi5 and its opposing activator, Cln3, are titrated against a constant number of sites on the genome occupied by the transcription factor SBF. This ‘titration of nuclear sites’ model captures all of the above data including the Whi5-independent increase in size of diploid cells. We also show that Whi5 dilution supports the START transition in this model but is not essential for proper cell size control. In summary, our modelling study suggests that the titration of nuclear sites is the dominant size-control mechanism in budding yeast. Because of the functional analogies between Cln3-Whi5-SBF in budding yeast and CycD-RB-E2F in mammalian cells, the titration of nuclear sites may be a conserved scenario in higher eukaryotes.

ABSTRACT. In proliferating cells, variability in cell size can rise from variability in growth rate, cell cycle length and the ratio between the two progeny cells. During recent years, the understanding about how single cell organisms restrict their size variability for the stable size distribution has built up rapidly. However, the mechanism of size control in proliferating mammalian cells remain largely unknown due to the lack of practical and accurate cell size measurement. Quantitative Phase Microscopy (QPM) is one of the gold standards to measure cell dry mass, whose sensitivity can be as low as 5 pg or 1% of the total cell dry mass. Here we turned this state-of-the-art technique into a robust and handy method to monitor the change of cell mass, with the throughput up to 1000 cells per minute and in different growth conditions. Using QPM, we observed the direct evidence of growth rate adjustment to reduce size variability in proliferating mammalian cells at two specific cell cycle stages. In late G1 and middle S phases, the correlation between cell mass and growth rate changes from positive to negative, thus the difference between large and small cells decreases. In the future, we will integrate the cell mass and growth rate quantification with other screenings to investigate the underlying mechanisms of growth rate dependent cell size control.

ABSTRACT. Robustness is a critical ability of biological oscillators to function in environmental perturbations. Although central architectures that support robust oscillations have been extensively studied, networks containing the same core vary drastically in their potential to oscillate, and it remains elusive what peripheral modifications to the core contribute to the large variation. We computationally generate a complete atlas of two- and three-node oscillators, to systematically analyze the association between network structures and robustness. We found that, while certain core topologies are essential for producing a robust oscillator, local structures can substantially modulate the degree of the robustness. Most strikingly, local nodes receiving incoherent (positive plus negative) or coherent (both positive or both negative) inputs promote or attenuate the overall network robustness significantly in an additive manner. These motifs are validated in larger-scale networks. Additionally, we found that incoherent inputs are enriched in almost all known natural and synthetic oscillators, suggesting that incoherent inputs may be a generalizable design principle that promotes oscillatory robustness. Our findings underscore the importance of local modifications besides robust cores, which explain why auxiliary structures not required for oscillation are evolutionarily conserved, and further suggest simple ways to evolve or design robust oscillators. Experimentally, we use microfluidics and fluorescence microscopy to investigate how network structures are linked to the essential functions of early embryonic cell cycles.

ABSTRACT. The brain has been studied as a complex networked system using various global and local network measures which can describe different aspects of the brain network. Among them, communicability is an extended measure of information flow by considering multiple paths between node pairs. The significance of communicability has been validated by many researchers studying the organizational principles of the brain, and also by clinicians dealing with lesion-like brain disorders such as stroke and multiple sclerosis. However, there is no study about the way of increasing communicability in damaged brain networks, which might provide us with a clue for brain rehabilitation. To tackle this problem, we have extracted structural brain networks of 40 normal adults from WU-Minn HCP's T1w and DTI image data. We then simulated brain disorders by attenuating or deleting some of the edges in the network. We further tested different edge addition strategies to find a method for restoring the communicability. As a result, we found that the optimal edge addition strategies depend on network attack methods. Among those, we found that there is a robust strategy for restoring the communicability regardless of the state of a damaged network. Our study provides a novel insight into the rehabilitation strategy for damaged brain networks in view of the network communicability.

ABSTRACT. Recently, the focus of systems neuroscience shifted to the determination of brain regions that allow the control of a whole-brain structural network. Finding the minimum dominating sets (MDSets) of regions, which potentially provide efficient sources of influence and information dispersal, can be a starting point for establishing a control strategy of the brain network. An intriguing question then arises as to whether such sets of regions have any important functional characteristics. In this study, we have investigated the MDSets of regions in the whole-brain structural network of human. In general, high-degree regions are more likely to be dominator nodes than low-degree regions, but, unexpectedly, the MDSets of regions include relatively low-degree regions and are minimally overlapped with rich-club regions. We further investigated the role of the MDSets of regions through network attack simulation and compared the result with that of the attack on the rich-club regions. We found that attacking the rich-club regions significantly decreases the global efficiency of the network while attacking the MDSets of regions significantly decreases the local efficiency of the network. Our study indicates that MDSets of regions might play a crucial role in local communication of the whole-brain structural network.

ABSTRACT. Cardiac action potentials (APs) are generated by the continuous dynamic balance between inward and outward voltage-gated ion channel currents. Gain or loss of channel function results in current imbalances, translating to AP prolongation (balance tipped in the inward) or shortening (balance tipped outward). Early afterdepolarizations (EADs), a known triggers for torsade de pointes arrhythmia, occur at a threshold increase in net inward current (1,2). Mild reduction in either or both of the outward IK1 and IKr currents due to mutations, results in higher susceptibility to EADs (3,4). We used the O’Hara-Rudy model of the undiseased human heart (5) to study the detailed mechanism by which EADs are generated in midmyocardial cells. As typical of dynamic instabilities, EADs result from a cascade of precipitating events. The AP transitions from a quasi- to a fully unstable state happens abruptly, for a particular value of current imbalance. Careful manipulation of the model was required to resolve these events, including complete abrogation of IK1, combined with mild hERG blockade. The results of our simulations suggest that EAD genesis is initiated by the spontaneous loss of IK1 at a critical inward current imbalance, due to prolongation of the AP at -50 mV (appearing as a shoulder in the waveform). IK1 plays a critical role in resetting the initial conditions of the system prior to the start of each AP cycle. Wherein the loss of this current results in: 1) abnormal buildup of Cav1.2 channels in the fast-inactivated state; 2) increase of the recovery-induced ICaL; and 3) a vicious circle of error propagation, consisting of abnormal Ca2+ release from internal stores, CAMKII activation, Cav1.2 phosphorylation, eventually increasing further Ca2+ influx. The accumulation of Ca2+, together with an increasing diastolic inward current generated by the Na+-Ca2+ exchanger, progressively moves the cell toward a more depolarized pro-arrhythmic state.

References

1. Weiss JN, Garfinkel A, Karagueuzian HS, Chen PS, Qu Z. Early afterdepolarizations and cardiac arrhythmias. Hear Rhythm. 2010; 2. Gilmour RF. Early afterdepolarization-induced triggered activity: Initiation and reinitiation of reentrant arrhythmias [Internet]. Vol. 1, Heart Rhythm. 2004 [cited 2017 Mar 29]. p. 449–50. Available from: http://www.sciencedirect.com/science/article/pii/S1547527104003923 3. Sartiani L, Bettiol E, Stillitano F, Mugelli A, Cerbai E, Jaconi ME. Developmental Changes in Cardiomyocytes Differentiated from Human Embryonic Stem Cells: A Molecular and Electrophysiological Approach. Stem Cells [Internet]. 2007 May [cited 2017 Mar 27];25(5):1136–44. Available from: http://www.ncbi.nlm.nih.gov/pubmed/17255522 4. Pogwizd SM, Schlotthauer K, Li L, Yuan W, Bers DM. Arrhythmogenesis and contractile dysfunction in heart failure: Roles of sodium-calcium exchange, inward rectifier potassium current, and residual beta-adrenergic responsiveness. Circ Res [Internet]. 2001 Jun 8 [cited 2017 Mar 28];88(11):1159–67. Available from: http://www.ncbi.nlm.nih.gov/pubmed/11397782 5. Rudy Y, Hara TO. Simulation of the Undiseased Human Cardiac Ventricular Action Potential : Model Formulation and Experimental Validation. 2011;7(5).

ABSTRACT. Functional Magnetic Resonance Imaging (fMRI) indirectly measures brain activity by registering changes in oxygen content that occur in response to neuronal activity. fMRI is being used in thousands of research articles and clinics each year, and is one of the two most common measurement techniques for recording brain activity. Despite this common usage, classical analysis of fMRI data is quite simplistic: it essentially consists of a correlation analysis with a pre-defined shape, even though it is known that the fMRI signal can behave very differently in different parts of the brain. This prompts development of new ways to interpret fMRI data that are more firmly based on our mechanistic understanding of how the signal is produced.

We have developed a mechanistic model based on ordinary differential equations that describes the intracellular signalling that connects neuronal glutamate and GABA release with vascular and metabolic mechanisms controlling brain oxygen level; this constitutes the previously unmodelled part of the basis for the fMRI signal. We have thereafter combined this model with the previously developed Balloon model, which describes the dynamic interplay between blood volume and flow. The combined model can describe fMRI data from several different clinical studies, featuring different types of response shapes in different parts of the brain, and can describe both estimation data used for parameter fitting, and independent validation data, used for model testing. Some of the key mechanisms in the model has also been tested using optogenetic experiments, where specific parts of the brain can be stimulated directly using light. Finally, using advanced uncertainty analysis, we have also characterized features in the model that can be uniquely identified from ordinary patient data, and that thus can serve as a new type of biologically based biomarkers, which can be used for patient stratification and diagnosis in a clinical setting.

ABSTRACT. The brain is a large complex network where various cognitive functions are embedded. To understand its dynamics, a system-level analysis is required. Resting-state recordings of the brain such as MEG or fMRI have been used as a common tool for studying resting-state brain networks and they show time-varying interactions across brain regions for specific cognitive functions. So far, Pearson’s correlation coefficient has been primarily used to measure the inter-regional relationships, but complex interactions among multiple brain regions cannot be fully quantified by this due to the assumption of the Pearson’s correlation that the pairwise regional interactions are independent of each other. Hence, in this study, we introduced the maximum entropy model based on resting-state MEG recordings to characterize the complex interdependent interactions between regions. Using this model, we further carried out energy landscape analysis to find out distinct resting-state dynamics for eight different frequency bands. We could observe state transitions upon each frequency-dependent landscape and found different dynamic activity patterns for each specific resting-state network. We also found local minima of each landscape and identified regions overlapped with the minima of multiple landscapes, which indicates control target regions for cross-frequency coupling. In summary, our study provides a novel insight into the brain resting-state dynamics from energy landscape analysis of frequency-dependent brain regional activities obtained by high temporal resolution of MEG.