- Email: [email protected]
Ensemble control for cell cycle synchronization of heterogeneous cell populations⁎
Proceedings, Proceedings, 7th 7th IFAC IFAC Conference Conference on on Proceedings, 7th IFAC Conference Available online at www.sciencedirect.com Foundations of Systems Biology Engineering Proceedings, 7th IFAC Conference on Foundations of Systems Biology in inon Engineering Proceedings, 7th IFAC Conference on Foundations of Systems Biology Engineering Chicago, Illinois, USA, 5-8, Foundations of Systems Biology in2018 Engineering Chicago, Illinois, USA, August August 5-8,in 2018 Foundations of Systems Biology Engineering Chicago, USA, 5-8, Chicago, Illinois, Illinois, USA, August August 5-8,in2018 2018 Chicago, Illinois, USA, August 5-8, 2018
ScienceDirect
IFAC PapersOnLine 51-19 (2018) 44–47
Ensemble control for cell cycle synchronization Ensemble Ensemble control control for for cell cell cycle cycle synchronization synchronization of heterogeneous heterogeneous cell populations Ensemble control for cell cycle synchronization of cell populations of heterogeneous cell populations ¨ of heterogeneous cell populations Karsten Kuritz Dyck Karsten Kuritz Dirke Dirke Imig Imig Michael Michael Dyck Frank Frank Allg Allgower ower ¨
Karsten ¨¨ Karsten Kuritz Kuritz Dirke Dirke Imig Imig Michael Michael Dyck Dyck Frank Frank Allg Allgower ower Karsten Kuritz Dirke Imig Michael Dyck Frank Allgower ¨ Institute for Systems Theory and Automatic Control, Institute for Systems Theory and Automatic Control, Institute for Systems Theory and Automatic Control, University of Stuttgart, Stuttgart, Pfaffenwaldring 9,Automatic 70569 Stuttgart, Stuttgart, Germany Institute for Systems Theory and9, Control,Germany University of Pfaffenwaldring 70569 Institute for Systems Theory and9, Control,Germany University of Pfaffenwaldring 70569 (e-mail: [email protected]) University of Stuttgart, Stuttgart, Pfaffenwaldring 9,Automatic 70569 Stuttgart, Stuttgart, Germany (e-mail: [email protected]) University of Stuttgart, 9, 70569 Stuttgart, Germany (e-mail: [email protected]) (e-mail: Pfaffenwaldring [email protected]) (e-mail: [email protected]) Abstract: We We investigate investigate the the problem problem of of synchronizing synchronizing aa population population of of cellular cellular oscillators oscillators in in Abstract: Abstract: We the of aa population of oscillators in their cell cell cycle. cycle. Restrictions onproblem the observability observability and controllability controllability of the the population imposed Abstract: We investigate investigate theon problem of synchronizing synchronizing population of cellular cellular oscillators in their Restrictions the and of population imposed Abstract: We investigate theon problem oftosynchronizing a population of cellular oscillators in their cell cycle. Restrictions the observability and controllability of the population imposed by the nature of cell biology give rise an ensemble control problem specified by finding their cell cycle. Restrictions on the observability and controllability of the population imposed by the nature of cell biology give rise to an ensemble control problem specified byimposed finding their cell cycle. Restrictions on the observability and controllability of the population by the nature of cell biology give rise to an ensemble control problem specified by finding a broadcast broadcast input based on the the distribution ofensemble the population. population. We recently showed, that the the nature of cell biology give rise to anof control We problem specified by that finding aby input based on distribution the recently showed, the bybroadcast the nature cell biology give rise to an control from problem specified by that finding a based the distribution of the population. We recently showed, the problem caninput beofsolved solved byon a passivity-based passivity-based control law derived the reduced reduced phase model a broadcast input based on the distribution ofensemble the law population. We recently showed, that the problem can be by a control derived from the phase model a broadcast based the distribution of the population. We recently showed, the problem can be solved by aa passivity-based control law derived from the reduced phase model representation of the population. In the the present present paper we investigate investigate the performance of the the problem caninput beof solved byon passivity-based control law derived from the reduced phasethat model representation the population. In paper we the performance of problem can be solved by a passivity-based control law derived from the reduced phase model representation the paper investigate performance the control law law in in of anthe in population. silico study study In where we simulate simulate heterogeneous cell populationof with representation of the population. In the present present paper aawe we investigate the thecell performance ofwith the control an in silico where we heterogeneous population representation population. In the present paperthe investigate performance ofwith the control law an in silico where we simulate aawe heterogeneous cell population an individual-based simulation framework. Therein, dynamics ofthe each cell are captured captured control law in in of anthe in simulation silico study studyframework. where we Therein, simulate heterogeneous cell population with an individual-based the dynamics of each cell are control law in an in silico study where we simulate a heterogeneous cell population with an individual-based simulation framework. Therein, the dynamics of each cell are captured by ordinary differential equations. Resampling of cellular states at cell division reflects the an individual-based simulation framework. Therein, the dynamics each cell are captured by ordinary differential equations. Resampling of cellular states at of cell division reflects the an individual-based simulation framework. Therein, the dynamics each are captured by ordinary differential equations. of cellular states at cell division reflects the naturally occurring heterogeneity in Resampling a population population due to uneven distribution ofcell proteins during by ordinary differential equations. Resampling of cellular states at of cell division reflects the naturally occurring heterogeneity in a due to uneven distribution of proteins during by ordinary differential equations. Resampling of cellular states atscenarios cell division reflects the naturally occurring heterogeneity in a population due to uneven distribution of proteins during cytokinesis. Performance of the controller is examined in different where we vary naturally occurring heterogeneity in a population due to uneven distribution of proteins during cytokinesis. Performance of the controller is examined in different scenarios where we vary naturally occurring heterogeneity in a population due to uneven distribution of proteins during cytokinesis. Performance of the controller is examined in different scenarios where we vary variance during during resampling against the input input strength. Our results indicate indicate suitable for cytokinesis. Performance ofagainst the controller is strength. examinedOur in different scenarios where range we vary variance resampling the results aa suitable for cytokinesis. Performance ofagainst the is controller isdespite examined in different scenarios where range we vary variance during resampling the strength. results indicate aaheterogeneity. suitable range for input strength where synchrony achieved theOur naturally occurring heterogeneity. variance during resampling against the input input strength. Our resultsoccurring indicate suitable range for input strength where synchrony is achieved despite the naturally variance duringwhere resampling against the inputdespite strength. resultsoccurring indicate aheterogeneity. suitable range for input strength synchrony is the naturally input strength where synchrony is achieved achieved despite theOur naturally occurring heterogeneity. input strength where synchrony is achieved despite the Hosting naturally occurring © 2018, IFAC (International Federation of Automatic Control) by Elsevier Ltd.heterogeneity. All rights reserved. Keywords: cell cell cycle, cycle, passivity, passivity, ensemble ensemble control, control, multi-agent multi-agent systems, systems, phase phase response response Keywords: Keywords: cell cycle, passivity, ensemble control, multi-agent systems, phase response Keywords: cell cycle, passivity, ensemble control, multi-agent systems, phase response Keywords: cell cycle, passivity, ensemble control, multi-agent systems, phase response 1. INTRODUCTION INTRODUCTION more realistic realistic experimental experimental observation observation is is composed composed of of 1. more 1. INTRODUCTION INTRODUCTION more realistic observation is of representative samples drawn from the population from 1. more realistic experimental experimental observation is composed composed of representative samples drawn from the population from The cell cell cycle cycle is is central central to life life and and the the mechanisms mechanisms which which more 1. INTRODUCTION realistic experimental observation iscycle composed of The to representative samples drawn from the population from which the distribution of cells in the cell must be representative samples drawn from the population from which the distribution of cells from in the cell cycle must be The cell cell cycle is central central to life lifecycles and the the mechanisms which give rise to possibly infinite of reproduction have The cycle is to and mechanisms which representative samples drawn the population from give rise to possibly infinite cycles of reproduction have which the distribution of cells in the cell cycle must be reconstructed. 2) Only Only of broadcast input signals can be which the distribution cells in the cellsignals cycle must be The cell cycle is central toare lifecycles and the mechanisms which 2) broadcast input can give rise to studied possibly infinite ofthe reproduction have reconstructed. long been and still focus give rise possibly infinite reproduction have which the distribution cells in the cellsignals cycle must be long beento studied and are cycles still in inof the focus of of ongoongoreconstructed. 2) broadcast input can realized. reconstructed. 2) Only Only of broadcast input signals can be be give rise to possibly infinite cycles of reproduction have realized. long been studied and are still in the focus of ongoing research. Besides its central role in reproduction, long been studied and are still in the focus of ongoreconstructed. 2) Only broadcast input signals can be ing research. Besides itsare central role infocus reproduction, We realized. longresearch. been studied andits still in the of ongo- realized. We recently recently proposed proposed aa passivity-based passivity-based control control law law for for ing research. Besides its central central role in reproduction, reproduction, tissue growth and various diseases including ing role in tissue growth Besides and renewal, renewal, various diseases including realized. We recently proposed a passivity-based control law for a population of identical cellular oscillators which We recently proposed a passivity-based control lawsynfor ing research. Besides its central role in reproduction, aWe population of identical cellular oscillators which syntissue growth and renewal, various diseases including Alzheimer’s disease and cancer are caused by malfunction tissue growth and renewal, various diseases including recently proposed a passivity-based control law for Alzheimer’s disease and cancer are caused by malfunction achronizes population of identical cellular oscillators which synchronizes the agents in their cell cycle (Kuritz et al., 2018). a population of identical cellular oscillators which syntissue growth and renewal, various diseases including the agents in their cell cycle (Kuritz et al., 2018). Alzheimer’s diseasecontrolled and cancer cancerprocess are caused caused by malfunction malfunction within this highly (Zhivotovsky and Alzheimer’s disease and are by achronizes population of identical cellular oscillators which synwithin this highly controlled process (Zhivotovsky and the agents in their cell cycle (Kuritz et al., 2018). The controller is based on a PDE of the reduced phase chronizes the agents in their cycle et al., phase 2018). Alzheimer’s disease and canceraiming are caused malfunction controller is based on acell PDE of (Kuritz the reduced within this2010). highly controlled process (Zhivotovsky and The Orrenius, Approaches at aa by deeper within this highly controlled process and chronizes the agents inof their cycle et1). al.,In 2018). Orrenius, 2010). Approaches aiming at(Zhivotovsky deeper underunderThe controller controller is based based on PDE of (Kuritz the (Fig. reduced phase model representation the population this The is on aacell PDE of the reduced phase within this highly controlled process (Zhivotovsky and model representation of the population (Fig. 1). In this Orrenius, 2010). Approaches aiming at a deeper understanding of the cell cycle machinery of controlOrrenius, 2010). Approaches aimingand at aways deeper under- The controller is based on a PDE of the reduced phase standing of the cell cycle machinery and ways of controlmodel representation of the population (Fig. 1). In this idealized description, the complex dynamics of oscillators model representation of the population (Fig. 1). In this Orrenius, 2010). Approaches aiming at a deeper underidealized description, the complex dynamics of oscillators standing of the cell cycle machinery and ways of controlling it are therefor crucial for improving the treatment of standing of the cell cycle machinery and ways of controlmodel representation of the population (Fig. 1). In this ling it are therefor crucial for improving the treatment of idealized description, the complex dynamics of oscillators on a limit cycle are captured by a single dynamic on the idealized description, the complex dynamics of oscillators standing of the cell cycle machinery and ways of controlon a limit cycle are captured by a single dynamic on the ling it it diseases. are therefor therefor crucial crucial for for improving improving the the treatment treatment of of idealized description, the complex dynamics of oscillators many ling are many diseases. on a limit cycle are captured by a single dynamic on the unit circle. The theory is based on the concept of weakly on a limit cycle are captured by a single dynamic on the ling it are therefor crucial for improving the treatment of unit circle. The theory is based on the concept of weakly many diseases. Mathematically, many diseases. the a circle. limitsystems cycle are(Hoppensteadt captured byon a and single dynamic on the Mathematically, the cell cell cycle cycle machinery machinery can can be be modeled modeled on unit The theory is based the concept of weakly coupled Izhikevich, 1997) unit circle. The theory is based on the concept of weakly many diseases. systems (Hoppensteadt and Izhikevich, 1997) Mathematically, the cell cellby cycle machinery can be beequations modeled coupled as aa dynamical system ordinary differential Mathematically, the cycle machinery can modeled unit circle. The theory is based on the concept of weakly as dynamical system by ordinary differential equations coupled systems (Hoppensteadt and Izhikevich, 1997) which imposes certain constraints input systems (Hoppensteadt andthe Izhikevich, 1997) Mathematically, the cell cycle machinery can beequations modeled coupled which imposes certain constraints on on the input strength. strength. as dynamical system by ordinary differential equations (ODEs) of form as aa dynamical system by ordinary differential coupled systems (Hoppensteadt and Izhikevich, 1997) (ODEs) of the the general general form which imposes certain constraints on the input strength. Synchronization of oscillators has been studied in differwhich imposes certain constraints on the input strength. as a dynamical system by ordinary differential equations Synchronization of oscillators has been studied in differ(ODEs) of of the the general general form form (ODEs) which imposes certain constraints on thestudied input strength. ˙ x = f (x, u) . (1) Synchronization of oscillators has been in different areas. Budding yeast synchronization via face locking ˙ Synchronization of oscillators has been studied in differx = f (x, u) . (1) (ODEs) of the general form areas. Budding yeast synchronization via faceinlocking Synchronization of oscillators has been studied differ= ff (x, (x, u) u) .. (1) ent xx˙˙ = (1) ent areas. Budding yeast synchronization via locking can be achieved by aa periodic extrinsic signal areas. Budding yeast synchronization via face face(Charvin locking Therein, molecular can be achieved by periodic extrinsic signal (Charvin Therein, the the states states xx represent represent different molecular species species x˙ = f (x,different u) . (1) ent ent areas. Budding yeast synchronization via face locking can be achieved by a periodic extrinsic signal (Charvin Therein, thewhich states xxcan represent different molecularinputs species et Similarly, diffusive coupling between oscilbe2009). achieved by a periodic (Charvin in the be by u Therein, the states represent different molecular species et al., al., 2009). Similarly, diffusiveextrinsic couplingsignal between oscilin the cell cell which can be affected affected by external external inputs u can becauses achieved bycertain a periodic extrinsic signal (Charvin Therein, thewhich states xcan represent different molecularinputs species et al., 2009). Similarly, diffusive coupling between oscilin theas cell which can be affected by external inputs u can lators under conditions synchronization of et al., 2009). Similarly, diffusive coupling between oscilsuch media, drugs, optogenetic approaches or enviin the cell be affected by external u such as media, drugs, optogenetic approaches or envilators causes under certain conditions synchronization of et al., 2009). Similarly, diffusive coupling between oscilin the cell which can be affected by external inputs u lators causes under certain conditions synchronization of such as media, drugs, optogenetic approaches or envia population (Montenbruck et al., 2015). While diffusive lators causes under certain conditions synchronization of ronmental factors (Milias-Argeitis et al., 2016). Besides such as media, drugs, optogenetic et approaches or envi- a population (Montenbruck et al., 2015). While diffusive ronmental factors (Milias-Argeitis al., 2016). Besides causes under certain conditions synchronization of such as media, drugs, optogenetic approaches or envi- lators aacoupling population (Montenbruck et al., 2015). While diffusive ronmental factors (Milias-Argeitis et al., 2016). Besides of cell cycle dynamics has not been observed, population (Montenbruck et al., 2015). While diffusive an agent-based description where each agent represents ronmental factors (Milias-Argeitis et al., 2016). Besides coupling of cell cycle dynamics has not While been observed, an agent-based description where each agent represents a population (Montenbruck et al., 2015). diffusive ronmental factors (Milias-Argeitis et al., 2016). Besides coupling of cell cycle dynamics has not been observed, an agent-based description where each agent represents open entrainment to oscillator is cell cycle dynamics has not been observed, aa single cell with dynamics Eq. an agent-based description each given agent by represents open loop loopof to an an external external oscillator is comcomsingle cell oscillating oscillating withwhere dynamics given by Eq. (1), (1), coupling ofentrainment cell cycle dynamics has been observed, ansingle agent-based description each agentwith represents open entrainment to external oscillator is coma single cell oscillating oscillating withwhere dynamics given by Eq. (1), coupling monlyloop observed. However, open loopnot approaches might open loop entrainment to an an external oscillator ismight comcell populations can be modeled partial aproliferating cell with dynamics given by Eq. (1), monly observed. However, open loop approaches proliferating cell populations can be modeled with partial open loop entrainment to an external oscillator is coma single cell oscillating with dynamics given by Eq. (1), monly observed. However, open loop loop approaches might proliferating cell populations can be modeled with partial result in anti-phase synchronization of cellular oscillators monly observed. However, open approaches might differential equations (PDEs) which describe the evoproliferating cell populations can be modeled with partial differential equations (PDEs) which describe the evoresult in anti-phase synchronization of cellular oscillators monly observed. However, open loop approaches might proliferating cell populations can be modeled with partial result in anti-phase of oscillators differential equations (PDEs) which describe the evoevowhile our proposed closed assures phase in anti-phase synchronization of cellular cellular oscillators lution of number density of The of differential the while our proposedsynchronization closed loop loop feedback feedback assures phase lution of the theequations number (PDEs) density which of cells. cells.describe The problem problem of result result in anti-phase synchronization of cellular oscillators differential equations (PDEs) which describe the evowhile our proposed closed loop feedback assures phase lution of the number density of cells. The problem of synchronization of the whole population. Here we present while our proposed closed loop feedback assures phase synchronizing a cell population within their cell cycle is lution of the number density of cells. The problem of synchronization of the whole population. Here we present synchronizing a cell population within their cell cycle is while our proposed closed loop feedback assures phase lution of the number density of cells. The problem of synchronization of the whole population. Here we present synchronizing a cell population within their cell cycle is a proof of concept study, where we apply the previously of the whole population. Here we present subject to several constraints imposed by the nature synchronizing a cell population within their cell cycle of is asynchronization proof of concept study, where we apply the previously subject to several constraints imposed by the nature of synchronization of the whole population. Here we present synchronizing a cell population within their cell cycle is proof study, where the previously subject to several several constraintsobservation imposed by byof the nature of aadeveloped ensemble control law to realistic proof of of concept concept study, where we apply the population previously cell biology. 1) the cell subject to constraints imposed of developed ensemble control lawwe to aaapply realistic cell biology. 1) Experimental Experimental observation ofthe the nature cell cycle cycle adeveloped proofheterogeneity of concept study, where the population previously subject to several constraints imposed byof the nature of developed ensemble control lawwe to by aapply realistic population cell biology. 1) Experimental observation of the cell cycle where is introduced resampling of ensemble control law to a realistic population state of individual cells over time is barely possible. A cell biology. 1) Experimental observation the cell cycle state of individual cells overobservation time is barely possible. A where heterogeneity is introduced by resampling of cell cell developed ensemble control law to a realistic population cell biology. 1) Experimental of the cell cycle where heterogeneity is introduced by resampling of state of individual cells over time is barely possible. A components at cell division as unequal partitioning of where heterogeneity is introduced by resampling of cell cell state of individual cells over time is barely possible. A components at cell division as unequal partitioning of The authors thank the German Research Foundation (DFG) for where heterogeneity is introduced by resampling of cell state individual cells over Research time is Foundation barely possible. A components The of authors thank the German (DFG) for at cell division as unequal partitioning of proteins after cell division plays a major role for cell-tocomponents at cell division as unequal partitioning of proteins after cell division plays a major role for cell-to The authors authors thank the German Research Foundation (DFG) and for financial supportthank of the thethe project under grant number number AL316/14-1 and German Research Foundation (DFG) for components atcell cell division as unequal partitioning of financial support of project under grant AL316/14-1 proteins after division plays a major role for cell-to The cell heterogeneity (Huh and Paulsson, 2011). The obtained proteins after cell division plays a major role for cell-toThe authors thank the German Research Foundation (DFG) for cell heterogeneity (Huh and Paulsson, 2011). The obtained financial support of the project under grant number AL316/14-1 and within the Cluster of Excellence in Simulation Technology (EXC 310/2) financial project in under grant number AL316/14-1 and proteins after cell division plays a major role for cell-towithin thesupport Clusterof of the Excellence Simulation Technology (EXC 310/2) cell heterogeneity (Huh and Paulsson, 2011). The obtained results indicate existence of strength cell heterogeneity (Huh and Paulsson, 2011).input The obtained financial the project in under grant number AL316/14-1 and within the Cluster of Excellence Simulation Technology (EXC results indicate the the existence of aa suitable suitable strength at the the University University ofof Stuttgart. within thesupport Clusterof of Excellence in Simulation Technology (EXC 310/2) 310/2) at Stuttgart. cell heterogeneity (Huh and Paulsson, 2011).input The obtained results indicate existence of input strength within the Clusterof Excellence in Simulation Technology (EXC 310/2) results indicate the the existence of aa suitable suitable input strength at the Stuttgart. at the University University ofof Stuttgart. results indicate the existence of a suitable input strength 2405-8963 © 2018, (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. at the University of IFAC Stuttgart.
Copyright 2018 44 Peer review© responsibility of International Federation of Automatic Copyright © under 2018 IFAC IFAC 44 Control. Copyright © 44 10.1016/j.ifacol.2018.09.034 Copyright © 2018 2018 IFAC IFAC 44 Copyright © 2018 IFAC 44
IFAC FOSBE 2018 Chicago, Illinois, USA, August 5-8, 2018
Karsten Kuritz et al. / IFAC PapersOnLine 51-19 (2018) 44–47
|m1 | = r ∈ [0, 1] is the order parameter which is equal to one if the population is synchronized (Kuramoto, 1984). The passivity-based controller that we recently developed provides a broadcast input signal based on a measure of the distribution of cells on the unit circle u = −Z, dp p( · , t) . (4) p We define d (θ) = d(p( · , t), θ) with ˆ ∗1 e− i θ + lnπ2 m ˆ 1m ˆ ∗1 − ln2π2 + 1 mˆ 1 ei θ , dp (θ) = − ln2π2 − 1 m
Mξ(t) x(t) ξ(t)
ψ = ψ1 ◦ ψ2
ξ(0) Mξ(0)
W γ
x(0)
θ(t)
θ(0)
ψ1
ψ2
S1
ψ1 1
and the distribution p(θ, t) in Eq. (4) is derived by transforming the distribution of cells obtained as measurement output such that it is smooth on the unit circle 2π −1 θ θ 2 2π ρ(θ, t) d θ . p(θ, t) = 2 2π ρ(θ, t)
Fig. 1. A neighborhood W of the limit cycle γ of an oscillator is invariantly foliated by isochrons Mξ . The flow maps isochrons to isochrons. The function ψ = ψ1 ◦ ψ2 maps an oscillator x(t) ∈ W uniquely to its phase on the unit circle θ(t) = ψ(x(t)).
0
window where synchronization is achieved. The paper is structured as follows: We first review the theoretical basis of the passivity-based ensemble control approach. We then describe the implementation of the controller in an individual-based simulation framework for heterogeneous populations. Subsequently, performance of the controller is studied in different scenarios where we varying input strength and noise during division.
A detailed derivation of the control law including proofs for convergence and controllability can be found in Kuritz et al. (2018). The aim of the present study is to evaluate the control law under realistic conditions where 1) cell division is subject to noise and 2) input acts on a state of the original system Eq. (1) rather than the reduced phase model Eq. (2). We are in particular interested whether there exists a biologically plausible setup in which the input can counteract division noise to achieve synchronization without interrupting the oscillations.
2. THEORETICAL BACKGROUND The notion of reduced phase models greatly simplifies the system to be controlled. Hence we review the basic concept of reduced phase models and phase response curves briefly (Hoppensteadt and Izhikevich, 1997; Kuritz et al., 2018). Given a system of the form of Eq. (1) with an exponentially stable limit cycle γ ∈ Rn with period T, then ˙ = ω + Z(θ) u , θ(t) ∈ S1 = [0, 2π) θ(t) (2)
3. IMPLEMENTATION We will first describe the dynamical model and the necessary steps to calculate the input from the full state model (Eq. (1)). We then outline the implementation of the controller in an individual-based simulation framework for heterogeneous populations (Imig et al., 2015). We used a five-variable skeleton model of the mammalian cell cycle by G´erard and Goldbeter (2011) consisting of the four main cyclin/Cdk complexes, the transcription factor E2F and the protein Cdc20 (Fig. 2). In the model, growth factors ensure the synthesis of the cyclin D/Cdk4–6 complex which rapidly reaches a steady state and promotes progression in the G1 phase by activating the transcription factor E2F. E2F causes activation of cyclin E/Cdk2 at the G1/S transition, and cyclin A/Cdk2 during S phase. During G2, cyclin A/Cdk2 triggers the activation of cyclin B/Cdk1, which leads to the G2/M transition. During mitosis, cyclin B/Cdk1 activates the protein Cdc20. Cdc20 creates a negative feedback loop involving cyclin A/Cdk2 and cyclin B/Cdk1 by promoting the degradation of these complexes. This negative feedback loop together with the positive loop of cyclin E/Cdk2 via E2F results in stable oscillations (Fig. 2 (b)). For the purpose of controlling the cell cycle progression, we extended the dynamics of the cyclin A/Cdk2 complex by an additive input d Ma Ma8 = fMA (x) + 1.6 1 − 8 u , (5) dt 4 + Ma8 which can be thought of e.g. an optogenetic signal causing a direct induction of cyclin A expression. At high levels of cyclin A the effect of the optogenetic input rapidly declines mimicking saturation of the amount of cyclin A. Strength of this input is tuned by the factor which is essential in order to keep the oscillators in a neighborhood W of the limit cycle γ. For the calculation of the input u(t) the phase of each agent has to be determined. This process follows
is a local canonical model for such an oscillator. The phase θ(t) of an oscillator proceeds with constant angular velocity ω = 2π T and an additional term where the input u acts on the phase via the phase response curve Z(θ). The phase response curve (PRC) describes the magnitude of phase changes after perturbing an oscillatory system. Reduced phase models and the concept of PRCs are based on the notion of isochrons introduced by Winfree (1974) (Fig. 1). Therein, the set of all initial conditions x(0) ∈ Rn of which the solution x(t) approaches the solution ξ(t), with ξ(0) ∈ γ is called an isochron of ξ(0). The flow maps isochrons to isochrons and Guckenheimer (1975) showed, that there always exists a neighborhood W of a limit cycle that is invariantly foliated by isochrons. The function ψ2 : W → γ maps a point in the neighborhood x ∈ Mξ ⊂ W to the generator of its isochron ξ ∈ γ. Furthermore, the periodic orbit of an oscillator is homeomorphic to the unit circle and the function ψ1 : γ → S1 maps the solution ξ(t) with ξ(0) ∈ γ to the solution of Eq. (2). The function ψ : W → S1 is a composition of ψ1 and ψ2 , mapping x(t) ∈ W uniquely to its corresponding phase θ(t) of the reduced phase model (Fig. 1). Synchrony of a population of oscillators can be measured by means of the first circular moment which is defined as 2π N ˆ 1 (t) = m1 (t) = ei θk (t) , or m ei θ ρ(θ, t) d θ , (3) k=1
45
0
for individual agents θk (t) or the distribution of agents on the unit circle ρ(θ, t), respectively. In the complex number m1 = mˆ1 = r ei φ , φ ∈ S1 is the mean phase and 45
IFAC FOSBE 2018 46 Chicago, Illinois, USA, August 5-8, 2018
Karsten Kuritz et al. / IFAC PapersOnLine 51-19 (2018) 44–47
• Upon cell division, heterogeneity is introduced by sampling the initial condition of the two daughter cells from a logarithmic normal distribution with sampling variance σ and the mean equal to the state values at the limit cycle γ at cell division. • The total number of cells is kept constant by removing a randomly determined cell from the list cells whenever a cell divides. A snapshot measurement is taken at the end of each simulation loop with a simulation time of tˆend = 1 h from which the input for the next loop is calculated.
(a)
(c)
3.0
cyclin B/Cdk1
2.5
G2/M ↖ 2.0
↙ M/G1
1.5 1.0
S/G2 ↗
0.5 0.0 0
1
2
cyclin A/Cdk2
3
d Ma dt d Mb dt d Me dt d E2F dt
Ma Kda + Ma Mb =vsb · Ma − Vdb · Cdc20 · Kdb + Mb Me =vse · E2F − Vde · Ma · Kde + Me E2Ftot − E2F =V1e2 f · K1e2 f + (E2Ftot − E2F)
=vsa · E2F − Vda · Cdc20 ·
d Cdc20 Cdc20tot − Cdc20 =V1dc20 · Mb · dt K1cdc20 + (Cdc20tot − Cdc20) − V2cdc20 ·
Cdc20 K2cdc20 + Cdc20
Fig. 2. (a) Scheme of the Cdk network driving the mammalian cell cycle (G´erard and Goldbeter, 2011). (b) Limit cycle oscillations shown as a projection onto the cyclin A/Cdk2 versus cyclin B/Cdk1 phase plane. (c) Corresponding system of ODEs, with complexes cyclin A/Cdk2, cyclin B/Cdk1 and cyclin E/Cdk2 abbreviated by Ma, Mb and Me, respectively.
Fig. 3. Scheme of the first layer of the simulation framework (adapted from Imig et al. (2015)).
the theory of reduced phase models, where ψ2 : W → γ maps a point x(t) ∈ W to the generator of its isochron ξ(t) ∈ γ (Fig. 1). However, due to the high complexity of the notion of isochrons we used a simplified mapping where ψˆ 2 maps the point x(t) to its closest point on the ˆ limit cycle ξ(t) ∈ γ. The relation ψ1 : γ → S1 which maps a point on the limit cycle γ to its corresponding position on the unit circle S1 which is proportional to cell age can easily be obtained by inversion of the simulated trajectory on the limit cycle ξ−1 (γ) = θ/ω or by employing ergodic principles as described in Kuritz et al. (2017). The composition ψ = ψ1 ◦ ψˆ 2 of the two mappings now maps the point x(t) ∈ W uniquely to its phase θ(t) ∈ S1 on the unit circle. The PRC Z was obtained by solving the appropriate adjoint equation using the dynamic modeling program XPPAUT (Malkin, 1956; Ermentrout, 2002). The distribution of cells on the unit circle was calculated by kernel density estimation from phase values θk using a von Mises distribution as kernel with manually tuned bandwidth. The input is calculated from the PRC and the distribution of cells on the unit circle by Eq. (4). The simulation framework is based on the individual-based population model by Imig et al. (2015) and was slightly adapted, resulting in the following framework (Fig. 3): • The simulation runs from time-points t0 until tend . • A list cells contains all initial mother-cells at t0 . • The cell with the shortest simulated time t0 is simulated until a specified time tˆend , or until the cell complies with a division condition.
4. SIMULATION RESULTS The cell population was initialized by applying the inverse : S1 → γ on samples drawn from a von mapping ψ−1 1 Mises distribution with circular moment m1 (t0 ) = 0.3. Simulations with varying input factor and sampling variance σ were run for 500 hours. Exemplary results obtained for sampling variance of 0.13 and 0.18 and input strength varying between 0 and 0.09 illustrate the different effects of the controller on the cell population (Fig. 4).
length first moment
(a)
σ = 0.13
1 0.8 0.6 0.4 0.2 0
(b)
length first moment
(b)
0
100 200 300 400 500
time in h
σ = 0.18
1 0.8 0.6 0.4 0.2 0
0
100 200 300 400 500
time in h
Fig. 4. Absolute value of the first circular moment for sampling variances σ = 0.13 (a) and σ = 0.18 (b). For < 0.02, the input is too small and the norm of the first circular moment remains constant or decreases, indicating a conversion towards a steady-state distribution. Increasing the input strength to = 0.05 results in stable 46
IFAC FOSBE 2018 Chicago, Illinois, USA, August 5-8, 2018
Karsten Kuritz et al. / IFAC PapersOnLine 51-19 (2018) 44–47
that the proposed control strategy is robust such that synchronization might also be achieved under stochastic dynamics. Typically synchronization of a cell population is achieved by environmental or chemical treatments. This procedure is often accompanied by complex side effects and the population loses synchrony rapidly. The presented controller can be used to culture cell populations which are synchronized while hardly affecting individual cells in their natural cell cycle. REFERENCES
oscillations with increasing length of the first moment. However, input strengths ≥ 0.09 are too strong causing the cells to diverge from their original limit cycle which results in loss of oscillation of parts of the population. Cell cycle arrest of a fraction of the population results in a very large amplitude of the length of the first moment. When the oscillating subpopulation passes by the arrested cells |m1 | increases. When the cells are on opposite sites of the cell cycle |m1 | is small. Such a behaviour can be seen for = 0.09 at the early time points with σ = 0.13 and throughout the whole simulation with σ = 0.18 (Fig. 4). Input strength and sampling variance were varied over a sufficiently large range, as illustrated in Fig. 5, to cover relevant outcomes and define regions where: 1) the input predominates over noise and the population synchronizes, 2) the population diverges and noise introduced during cell division predominates, or 3) the input is too strong causing loss of oscillations. The input strength which is necessary to synchronize the population increases with the variance until it reaches a point where the input causes loss of oscillations. The border of this region expands at higher variances to lower input strengths. A possible explanation for this observation is that higher variances during resampling result in larger deviations from the limit cycle during cell division such that lower inputs suffice to interrupt stable oscillations. The decline rate of |m1 | in the simulation with zero input and sampling variances of σ = 0.2 is comparable to a cell cycle length distribution with mean 24 h and variance 2 h. Such a distribution is in good accordance with the experimentally observed cell cycle lengths of mammalian cell cultures.
Charvin, G., Cross, F.R., and Siggia, E.D. (2009). Forced periodic expression of G1 cyclins phase-locks the budding yeast cell cycle. Proceedings of the National Academy of Sciences of the United States of America, 106, 6632–6637. doi:10.1073/pnas.0809227106. Ermentrout, B. (2002). Simulating, Analyzing, and Animating Dynamical Systems. Society for Industrial and Applied Mathematics. G´erard, C. and Goldbeter, A. (2011). A skeleton model for the network of cyclin-dependent kinases driving the mammalian cell cycle. Interface Focus, 1(1), 24–35. Guckenheimer, J. (1975). Isochrons and phaseless sets. J. Math. Biol., 1(3), 259–273. Hanahan, D. and Weinberg, R.A. (2011). Hallmarks of cancer: the next generation. Cell, 144(5), 646–74. Hoppensteadt, F.C. and Izhikevich, E.M. (1997). Weakly Connected Neural Networks, volume 126 of Applied Mathematical Sciences. Springer New York. Huh, D. and Paulsson, J. (2011). Non-genetic heterogeneity from stochastic partitioning at cell division. Nature Genetics, 43(2), 95–100. doi:10.1038/ng.729. Imig, D., Pollak, N., Strecker, T., Scheurich, P., Allgower, ¨ F., and Waldherr, S. (2015). An individual-based simulation framework for dynamic, heterogeneous cell populations during extrinsic stimulations. J. Coupled Syst. Multiscale Dyn., 3(2), 143–155. Kuramoto, Y. (1984). Chemical Oscillations, Waves, and Turbulence, volume 19 of Springer Series in Synergetics. Springer Berlin Heidelberg. Kuritz, K., Stohr, D., Pollak, N., and Allgower, F. (2017). ¨ ¨ On the relationship between cell cycle analysis with ergodic principles and age-structured cell population models. J. Theor. Biol., 414(7), 91–102. Kuritz, K., Halter, W., and Allgower, F. (2018). Passivity¨ based ensemble control for cell cycle synchronization. In R. Tempo, S. Yurkovich, and P. Misra (eds.), Emerg. Appl. Control Syst. Theory. Springer International Publishing, 1 edition. Malkin, I.G. (1956). Some Problems in Nonlinear Oscillation Theory. Gostexizdat, Moscow. Milias-Argeitis, A., Rullan, M., Aoki, S.K., Buchmann, P., and Khammash, M. (2016). Automated optogenetic feedback control for precise and robust regulation of gene expression and cell growth. Nat. Commun., 7, 12546. Montenbruck, J.M., Burger, M., and Allgower, F. (2015). ¨ ¨ Practical synchronization with diffusive couplings. Automatica, 53, 235–243. Winfree, A.T. (1974). Patterns of phase compromise in biological cycles. J. Math. Biol., 1(1), 73–93. Zhivotovsky, B. and Orrenius, S. (2010). Cell cycle and cell death in disease: past, present and future. J. Intern. Med., 268(5), 395–409.
oscilla
division noise predominates
tions
0.2 0.15 population synchronizes 0.1
0
0.02
0.04 0.06 input strength
0.08
upted
interr
resampling variance σ
0.3 0.25
47
0.1
Fig. 5. Parameter regions where we observe no synchronization due to predominant division noise, successful synchronization or loss of oscillations. 5. CONCLUSION In this paper we describe a proof-of-concept for the recently derived ensemble control approach to synchronize a cell population within the cell cycle (Kuritz et al., 2018). We expand our previous work which was derived for an idealized behavior of the cell population to more realistic conditions. To this end we used individual-based simulation framework where heterogeneity is introduced by resampling of states upon cell division (Imig et al., 2015). The simulation results indicate the existence of a suitable range in which the input strength can be varied while still achieving synchrony without interrupting the oscillations. Heterogeneity in the population can also be introduced by stochastic dynamics. We hypothesize 47
ScienceDirect
IFAC PapersOnLine 51-19 (2018) 44–47
Ensemble control for cell cycle synchronization Ensemble Ensemble control control for for cell cell cycle cycle synchronization synchronization of heterogeneous heterogeneous cell populations Ensemble control for cell cycle synchronization of cell populations of heterogeneous cell populations ¨ of heterogeneous cell populations Karsten Kuritz Dyck Karsten Kuritz Dirke Dirke Imig Imig Michael Michael Dyck Frank Frank Allg Allgower ower ¨
Karsten ¨¨ Karsten Kuritz Kuritz Dirke Dirke Imig Imig Michael Michael Dyck Dyck Frank Frank Allg Allgower ower Karsten Kuritz Dirke Imig Michael Dyck Frank Allgower ¨ Institute for Systems Theory and Automatic Control, Institute for Systems Theory and Automatic Control, Institute for Systems Theory and Automatic Control, University of Stuttgart, Stuttgart, Pfaffenwaldring 9,Automatic 70569 Stuttgart, Stuttgart, Germany Institute for Systems Theory and9, Control,Germany University of Pfaffenwaldring 70569 Institute for Systems Theory and9, Control,Germany University of Pfaffenwaldring 70569 (e-mail: [email protected]) University of Stuttgart, Stuttgart, Pfaffenwaldring 9,Automatic 70569 Stuttgart, Stuttgart, Germany (e-mail: [email protected]) University of Stuttgart, 9, 70569 Stuttgart, Germany (e-mail: [email protected]) (e-mail: Pfaffenwaldring [email protected]) (e-mail: [email protected]) Abstract: We We investigate investigate the the problem problem of of synchronizing synchronizing aa population population of of cellular cellular oscillators oscillators in in Abstract: Abstract: We the of aa population of oscillators in their cell cell cycle. cycle. Restrictions onproblem the observability observability and controllability controllability of the the population imposed Abstract: We investigate investigate theon problem of synchronizing synchronizing population of cellular cellular oscillators in their Restrictions the and of population imposed Abstract: We investigate theon problem oftosynchronizing a population of cellular oscillators in their cell cycle. Restrictions the observability and controllability of the population imposed by the nature of cell biology give rise an ensemble control problem specified by finding their cell cycle. Restrictions on the observability and controllability of the population imposed by the nature of cell biology give rise to an ensemble control problem specified byimposed finding their cell cycle. Restrictions on the observability and controllability of the population by the nature of cell biology give rise to an ensemble control problem specified by finding a broadcast broadcast input based on the the distribution ofensemble the population. population. We recently showed, that the the nature of cell biology give rise to anof control We problem specified by that finding aby input based on distribution the recently showed, the bybroadcast the nature cell biology give rise to an control from problem specified by that finding a based the distribution of the population. We recently showed, the problem caninput beofsolved solved byon a passivity-based passivity-based control law derived the reduced reduced phase model a broadcast input based on the distribution ofensemble the law population. We recently showed, that the problem can be by a control derived from the phase model a broadcast based the distribution of the population. We recently showed, the problem can be solved by aa passivity-based control law derived from the reduced phase model representation of the population. In the the present present paper we investigate investigate the performance of the the problem caninput beof solved byon passivity-based control law derived from the reduced phasethat model representation the population. In paper we the performance of problem can be solved by a passivity-based control law derived from the reduced phase model representation the paper investigate performance the control law law in in of anthe in population. silico study study In where we simulate simulate heterogeneous cell populationof with representation of the population. In the present present paper aawe we investigate the thecell performance ofwith the control an in silico where we heterogeneous population representation population. In the present paperthe investigate performance ofwith the control law an in silico where we simulate aawe heterogeneous cell population an individual-based simulation framework. Therein, dynamics ofthe each cell are captured captured control law in in of anthe in simulation silico study studyframework. where we Therein, simulate heterogeneous cell population with an individual-based the dynamics of each cell are control law in an in silico study where we simulate a heterogeneous cell population with an individual-based simulation framework. Therein, the dynamics of each cell are captured by ordinary differential equations. Resampling of cellular states at cell division reflects the an individual-based simulation framework. Therein, the dynamics each cell are captured by ordinary differential equations. Resampling of cellular states at of cell division reflects the an individual-based simulation framework. Therein, the dynamics each are captured by ordinary differential equations. of cellular states at cell division reflects the naturally occurring heterogeneity in Resampling a population population due to uneven distribution ofcell proteins during by ordinary differential equations. Resampling of cellular states at of cell division reflects the naturally occurring heterogeneity in a due to uneven distribution of proteins during by ordinary differential equations. Resampling of cellular states atscenarios cell division reflects the naturally occurring heterogeneity in a population due to uneven distribution of proteins during cytokinesis. Performance of the controller is examined in different where we vary naturally occurring heterogeneity in a population due to uneven distribution of proteins during cytokinesis. Performance of the controller is examined in different scenarios where we vary naturally occurring heterogeneity in a population due to uneven distribution of proteins during cytokinesis. Performance of the controller is examined in different scenarios where we vary variance during during resampling against the input input strength. Our results indicate indicate suitable for cytokinesis. Performance ofagainst the controller is strength. examinedOur in different scenarios where range we vary variance resampling the results aa suitable for cytokinesis. Performance ofagainst the is controller isdespite examined in different scenarios where range we vary variance during resampling the strength. results indicate aaheterogeneity. suitable range for input strength where synchrony achieved theOur naturally occurring heterogeneity. variance during resampling against the input input strength. Our resultsoccurring indicate suitable range for input strength where synchrony is achieved despite the naturally variance duringwhere resampling against the inputdespite strength. resultsoccurring indicate aheterogeneity. suitable range for input strength synchrony is the naturally input strength where synchrony is achieved achieved despite theOur naturally occurring heterogeneity. input strength where synchrony is achieved despite the Hosting naturally occurring © 2018, IFAC (International Federation of Automatic Control) by Elsevier Ltd.heterogeneity. All rights reserved. Keywords: cell cell cycle, cycle, passivity, passivity, ensemble ensemble control, control, multi-agent multi-agent systems, systems, phase phase response response Keywords: Keywords: cell cycle, passivity, ensemble control, multi-agent systems, phase response Keywords: cell cycle, passivity, ensemble control, multi-agent systems, phase response Keywords: cell cycle, passivity, ensemble control, multi-agent systems, phase response 1. INTRODUCTION INTRODUCTION more realistic realistic experimental experimental observation observation is is composed composed of of 1. more 1. INTRODUCTION INTRODUCTION more realistic observation is of representative samples drawn from the population from 1. more realistic experimental experimental observation is composed composed of representative samples drawn from the population from The cell cell cycle cycle is is central central to life life and and the the mechanisms mechanisms which which more 1. INTRODUCTION realistic experimental observation iscycle composed of The to representative samples drawn from the population from which the distribution of cells in the cell must be representative samples drawn from the population from which the distribution of cells from in the cell cycle must be The cell cell cycle is central central to life lifecycles and the the mechanisms which give rise to possibly infinite of reproduction have The cycle is to and mechanisms which representative samples drawn the population from give rise to possibly infinite cycles of reproduction have which the distribution of cells in the cell cycle must be reconstructed. 2) Only Only of broadcast input signals can be which the distribution cells in the cellsignals cycle must be The cell cycle is central toare lifecycles and the mechanisms which 2) broadcast input can give rise to studied possibly infinite ofthe reproduction have reconstructed. long been and still focus give rise possibly infinite reproduction have which the distribution cells in the cellsignals cycle must be long beento studied and are cycles still in inof the focus of of ongoongoreconstructed. 2) broadcast input can realized. reconstructed. 2) Only Only of broadcast input signals can be be give rise to possibly infinite cycles of reproduction have realized. long been studied and are still in the focus of ongoing research. Besides its central role in reproduction, long been studied and are still in the focus of ongoreconstructed. 2) Only broadcast input signals can be ing research. Besides itsare central role infocus reproduction, We realized. longresearch. been studied andits still in the of ongo- realized. We recently recently proposed proposed aa passivity-based passivity-based control control law law for for ing research. Besides its central central role in reproduction, reproduction, tissue growth and various diseases including ing role in tissue growth Besides and renewal, renewal, various diseases including realized. We recently proposed a passivity-based control law for a population of identical cellular oscillators which We recently proposed a passivity-based control lawsynfor ing research. Besides its central role in reproduction, aWe population of identical cellular oscillators which syntissue growth and renewal, various diseases including Alzheimer’s disease and cancer are caused by malfunction tissue growth and renewal, various diseases including recently proposed a passivity-based control law for Alzheimer’s disease and cancer are caused by malfunction achronizes population of identical cellular oscillators which synchronizes the agents in their cell cycle (Kuritz et al., 2018). a population of identical cellular oscillators which syntissue growth and renewal, various diseases including the agents in their cell cycle (Kuritz et al., 2018). Alzheimer’s diseasecontrolled and cancer cancerprocess are caused caused by malfunction malfunction within this highly (Zhivotovsky and Alzheimer’s disease and are by achronizes population of identical cellular oscillators which synwithin this highly controlled process (Zhivotovsky and the agents in their cell cycle (Kuritz et al., 2018). The controller is based on a PDE of the reduced phase chronizes the agents in their cycle et al., phase 2018). Alzheimer’s disease and canceraiming are caused malfunction controller is based on acell PDE of (Kuritz the reduced within this2010). highly controlled process (Zhivotovsky and The Orrenius, Approaches at aa by deeper within this highly controlled process and chronizes the agents inof their cycle et1). al.,In 2018). Orrenius, 2010). Approaches aiming at(Zhivotovsky deeper underunderThe controller controller is based based on PDE of (Kuritz the (Fig. reduced phase model representation the population this The is on aacell PDE of the reduced phase within this highly controlled process (Zhivotovsky and model representation of the population (Fig. 1). In this Orrenius, 2010). Approaches aiming at a deeper understanding of the cell cycle machinery of controlOrrenius, 2010). Approaches aimingand at aways deeper under- The controller is based on a PDE of the reduced phase standing of the cell cycle machinery and ways of controlmodel representation of the population (Fig. 1). In this idealized description, the complex dynamics of oscillators model representation of the population (Fig. 1). In this Orrenius, 2010). Approaches aiming at a deeper underidealized description, the complex dynamics of oscillators standing of the cell cycle machinery and ways of controlling it are therefor crucial for improving the treatment of standing of the cell cycle machinery and ways of controlmodel representation of the population (Fig. 1). In this ling it are therefor crucial for improving the treatment of idealized description, the complex dynamics of oscillators on a limit cycle are captured by a single dynamic on the idealized description, the complex dynamics of oscillators standing of the cell cycle machinery and ways of controlon a limit cycle are captured by a single dynamic on the ling it it diseases. are therefor therefor crucial crucial for for improving improving the the treatment treatment of of idealized description, the complex dynamics of oscillators many ling are many diseases. on a limit cycle are captured by a single dynamic on the unit circle. The theory is based on the concept of weakly on a limit cycle are captured by a single dynamic on the ling it are therefor crucial for improving the treatment of unit circle. The theory is based on the concept of weakly many diseases. Mathematically, many diseases. the a circle. limitsystems cycle are(Hoppensteadt captured byon a and single dynamic on the Mathematically, the cell cell cycle cycle machinery machinery can can be be modeled modeled on unit The theory is based the concept of weakly coupled Izhikevich, 1997) unit circle. The theory is based on the concept of weakly many diseases. systems (Hoppensteadt and Izhikevich, 1997) Mathematically, the cell cellby cycle machinery can be beequations modeled coupled as aa dynamical system ordinary differential Mathematically, the cycle machinery can modeled unit circle. The theory is based on the concept of weakly as dynamical system by ordinary differential equations coupled systems (Hoppensteadt and Izhikevich, 1997) which imposes certain constraints input systems (Hoppensteadt andthe Izhikevich, 1997) Mathematically, the cell cycle machinery can beequations modeled coupled which imposes certain constraints on on the input strength. strength. as dynamical system by ordinary differential equations (ODEs) of form as aa dynamical system by ordinary differential coupled systems (Hoppensteadt and Izhikevich, 1997) (ODEs) of the the general general form which imposes certain constraints on the input strength. Synchronization of oscillators has been studied in differwhich imposes certain constraints on the input strength. as a dynamical system by ordinary differential equations Synchronization of oscillators has been studied in differ(ODEs) of of the the general general form form (ODEs) which imposes certain constraints on thestudied input strength. ˙ x = f (x, u) . (1) Synchronization of oscillators has been in different areas. Budding yeast synchronization via face locking ˙ Synchronization of oscillators has been studied in differx = f (x, u) . (1) (ODEs) of the general form areas. Budding yeast synchronization via faceinlocking Synchronization of oscillators has been studied differ= ff (x, (x, u) u) .. (1) ent xx˙˙ = (1) ent areas. Budding yeast synchronization via locking can be achieved by aa periodic extrinsic signal areas. Budding yeast synchronization via face face(Charvin locking Therein, molecular can be achieved by periodic extrinsic signal (Charvin Therein, the the states states xx represent represent different molecular species species x˙ = f (x,different u) . (1) ent ent areas. Budding yeast synchronization via face locking can be achieved by a periodic extrinsic signal (Charvin Therein, thewhich states xxcan represent different molecularinputs species et Similarly, diffusive coupling between oscilbe2009). achieved by a periodic (Charvin in the be by u Therein, the states represent different molecular species et al., al., 2009). Similarly, diffusiveextrinsic couplingsignal between oscilin the cell cell which can be affected affected by external external inputs u can becauses achieved bycertain a periodic extrinsic signal (Charvin Therein, thewhich states xcan represent different molecularinputs species et al., 2009). Similarly, diffusive coupling between oscilin theas cell which can be affected by external inputs u can lators under conditions synchronization of et al., 2009). Similarly, diffusive coupling between oscilsuch media, drugs, optogenetic approaches or enviin the cell be affected by external u such as media, drugs, optogenetic approaches or envilators causes under certain conditions synchronization of et al., 2009). Similarly, diffusive coupling between oscilin the cell which can be affected by external inputs u lators causes under certain conditions synchronization of such as media, drugs, optogenetic approaches or envia population (Montenbruck et al., 2015). While diffusive lators causes under certain conditions synchronization of ronmental factors (Milias-Argeitis et al., 2016). Besides such as media, drugs, optogenetic et approaches or envi- a population (Montenbruck et al., 2015). While diffusive ronmental factors (Milias-Argeitis al., 2016). Besides causes under certain conditions synchronization of such as media, drugs, optogenetic approaches or envi- lators aacoupling population (Montenbruck et al., 2015). While diffusive ronmental factors (Milias-Argeitis et al., 2016). Besides of cell cycle dynamics has not been observed, population (Montenbruck et al., 2015). While diffusive an agent-based description where each agent represents ronmental factors (Milias-Argeitis et al., 2016). Besides coupling of cell cycle dynamics has not While been observed, an agent-based description where each agent represents a population (Montenbruck et al., 2015). diffusive ronmental factors (Milias-Argeitis et al., 2016). Besides coupling of cell cycle dynamics has not been observed, an agent-based description where each agent represents open entrainment to oscillator is cell cycle dynamics has not been observed, aa single cell with dynamics Eq. an agent-based description each given agent by represents open loop loopof to an an external external oscillator is comcomsingle cell oscillating oscillating withwhere dynamics given by Eq. (1), (1), coupling ofentrainment cell cycle dynamics has been observed, ansingle agent-based description each agentwith represents open entrainment to external oscillator is coma single cell oscillating oscillating withwhere dynamics given by Eq. (1), coupling monlyloop observed. However, open loopnot approaches might open loop entrainment to an an external oscillator ismight comcell populations can be modeled partial aproliferating cell with dynamics given by Eq. (1), monly observed. However, open loop approaches proliferating cell populations can be modeled with partial open loop entrainment to an external oscillator is coma single cell oscillating with dynamics given by Eq. (1), monly observed. However, open loop loop approaches might proliferating cell populations can be modeled with partial result in anti-phase synchronization of cellular oscillators monly observed. However, open approaches might differential equations (PDEs) which describe the evoproliferating cell populations can be modeled with partial differential equations (PDEs) which describe the evoresult in anti-phase synchronization of cellular oscillators monly observed. However, open loop approaches might proliferating cell populations can be modeled with partial result in anti-phase of oscillators differential equations (PDEs) which describe the evoevowhile our proposed closed assures phase in anti-phase synchronization of cellular cellular oscillators lution of number density of The of differential the while our proposedsynchronization closed loop loop feedback feedback assures phase lution of the theequations number (PDEs) density which of cells. cells.describe The problem problem of result result in anti-phase synchronization of cellular oscillators differential equations (PDEs) which describe the evowhile our proposed closed loop feedback assures phase lution of the number density of cells. The problem of synchronization of the whole population. Here we present while our proposed closed loop feedback assures phase synchronizing a cell population within their cell cycle is lution of the number density of cells. The problem of synchronization of the whole population. Here we present synchronizing a cell population within their cell cycle is while our proposed closed loop feedback assures phase lution of the number density of cells. The problem of synchronization of the whole population. Here we present synchronizing a cell population within their cell cycle is a proof of concept study, where we apply the previously of the whole population. Here we present subject to several constraints imposed by the nature synchronizing a cell population within their cell cycle of is asynchronization proof of concept study, where we apply the previously subject to several constraints imposed by the nature of synchronization of the whole population. Here we present synchronizing a cell population within their cell cycle is proof study, where the previously subject to several several constraintsobservation imposed by byof the nature of aadeveloped ensemble control law to realistic proof of of concept concept study, where we apply the population previously cell biology. 1) the cell subject to constraints imposed of developed ensemble control lawwe to aaapply realistic cell biology. 1) Experimental Experimental observation ofthe the nature cell cycle cycle adeveloped proofheterogeneity of concept study, where the population previously subject to several constraints imposed byof the nature of developed ensemble control lawwe to by aapply realistic population cell biology. 1) Experimental observation of the cell cycle where is introduced resampling of ensemble control law to a realistic population state of individual cells over time is barely possible. A cell biology. 1) Experimental observation the cell cycle state of individual cells overobservation time is barely possible. A where heterogeneity is introduced by resampling of cell cell developed ensemble control law to a realistic population cell biology. 1) Experimental of the cell cycle where heterogeneity is introduced by resampling of state of individual cells over time is barely possible. A components at cell division as unequal partitioning of where heterogeneity is introduced by resampling of cell cell state of individual cells over time is barely possible. A components at cell division as unequal partitioning of The authors thank the German Research Foundation (DFG) for where heterogeneity is introduced by resampling of cell state individual cells over Research time is Foundation barely possible. A components The of authors thank the German (DFG) for at cell division as unequal partitioning of proteins after cell division plays a major role for cell-tocomponents at cell division as unequal partitioning of proteins after cell division plays a major role for cell-to The authors authors thank the German Research Foundation (DFG) and for financial supportthank of the thethe project under grant number number AL316/14-1 and German Research Foundation (DFG) for components atcell cell division as unequal partitioning of financial support of project under grant AL316/14-1 proteins after division plays a major role for cell-to The cell heterogeneity (Huh and Paulsson, 2011). The obtained proteins after cell division plays a major role for cell-toThe authors thank the German Research Foundation (DFG) for cell heterogeneity (Huh and Paulsson, 2011). The obtained financial support of the project under grant number AL316/14-1 and within the Cluster of Excellence in Simulation Technology (EXC 310/2) financial project in under grant number AL316/14-1 and proteins after cell division plays a major role for cell-towithin thesupport Clusterof of the Excellence Simulation Technology (EXC 310/2) cell heterogeneity (Huh and Paulsson, 2011). The obtained results indicate existence of strength cell heterogeneity (Huh and Paulsson, 2011).input The obtained financial the project in under grant number AL316/14-1 and within the Cluster of Excellence Simulation Technology (EXC results indicate the the existence of aa suitable suitable strength at the the University University ofof Stuttgart. within thesupport Clusterof of Excellence in Simulation Technology (EXC 310/2) 310/2) at Stuttgart. cell heterogeneity (Huh and Paulsson, 2011).input The obtained results indicate existence of input strength within the Clusterof Excellence in Simulation Technology (EXC 310/2) results indicate the the existence of aa suitable suitable input strength at the Stuttgart. at the University University ofof Stuttgart. results indicate the existence of a suitable input strength 2405-8963 © 2018, (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. at the University of IFAC Stuttgart.
Copyright 2018 44 Peer review© responsibility of International Federation of Automatic Copyright © under 2018 IFAC IFAC 44 Control. Copyright © 44 10.1016/j.ifacol.2018.09.034 Copyright © 2018 2018 IFAC IFAC 44 Copyright © 2018 IFAC 44
IFAC FOSBE 2018 Chicago, Illinois, USA, August 5-8, 2018
Karsten Kuritz et al. / IFAC PapersOnLine 51-19 (2018) 44–47
|m1 | = r ∈ [0, 1] is the order parameter which is equal to one if the population is synchronized (Kuramoto, 1984). The passivity-based controller that we recently developed provides a broadcast input signal based on a measure of the distribution of cells on the unit circle u = −Z, dp p( · , t) . (4) p We define d (θ) = d(p( · , t), θ) with ˆ ∗1 e− i θ + lnπ2 m ˆ 1m ˆ ∗1 − ln2π2 + 1 mˆ 1 ei θ , dp (θ) = − ln2π2 − 1 m
Mξ(t) x(t) ξ(t)
ψ = ψ1 ◦ ψ2
ξ(0) Mξ(0)
W γ
x(0)
θ(t)
θ(0)
ψ1
ψ2
S1
ψ1 1
and the distribution p(θ, t) in Eq. (4) is derived by transforming the distribution of cells obtained as measurement output such that it is smooth on the unit circle 2π −1 θ θ 2 2π ρ(θ, t) d θ . p(θ, t) = 2 2π ρ(θ, t)
Fig. 1. A neighborhood W of the limit cycle γ of an oscillator is invariantly foliated by isochrons Mξ . The flow maps isochrons to isochrons. The function ψ = ψ1 ◦ ψ2 maps an oscillator x(t) ∈ W uniquely to its phase on the unit circle θ(t) = ψ(x(t)).
0
window where synchronization is achieved. The paper is structured as follows: We first review the theoretical basis of the passivity-based ensemble control approach. We then describe the implementation of the controller in an individual-based simulation framework for heterogeneous populations. Subsequently, performance of the controller is studied in different scenarios where we varying input strength and noise during division.
A detailed derivation of the control law including proofs for convergence and controllability can be found in Kuritz et al. (2018). The aim of the present study is to evaluate the control law under realistic conditions where 1) cell division is subject to noise and 2) input acts on a state of the original system Eq. (1) rather than the reduced phase model Eq. (2). We are in particular interested whether there exists a biologically plausible setup in which the input can counteract division noise to achieve synchronization without interrupting the oscillations.
2. THEORETICAL BACKGROUND The notion of reduced phase models greatly simplifies the system to be controlled. Hence we review the basic concept of reduced phase models and phase response curves briefly (Hoppensteadt and Izhikevich, 1997; Kuritz et al., 2018). Given a system of the form of Eq. (1) with an exponentially stable limit cycle γ ∈ Rn with period T, then ˙ = ω + Z(θ) u , θ(t) ∈ S1 = [0, 2π) θ(t) (2)
3. IMPLEMENTATION We will first describe the dynamical model and the necessary steps to calculate the input from the full state model (Eq. (1)). We then outline the implementation of the controller in an individual-based simulation framework for heterogeneous populations (Imig et al., 2015). We used a five-variable skeleton model of the mammalian cell cycle by G´erard and Goldbeter (2011) consisting of the four main cyclin/Cdk complexes, the transcription factor E2F and the protein Cdc20 (Fig. 2). In the model, growth factors ensure the synthesis of the cyclin D/Cdk4–6 complex which rapidly reaches a steady state and promotes progression in the G1 phase by activating the transcription factor E2F. E2F causes activation of cyclin E/Cdk2 at the G1/S transition, and cyclin A/Cdk2 during S phase. During G2, cyclin A/Cdk2 triggers the activation of cyclin B/Cdk1, which leads to the G2/M transition. During mitosis, cyclin B/Cdk1 activates the protein Cdc20. Cdc20 creates a negative feedback loop involving cyclin A/Cdk2 and cyclin B/Cdk1 by promoting the degradation of these complexes. This negative feedback loop together with the positive loop of cyclin E/Cdk2 via E2F results in stable oscillations (Fig. 2 (b)). For the purpose of controlling the cell cycle progression, we extended the dynamics of the cyclin A/Cdk2 complex by an additive input d Ma Ma8 = fMA (x) + 1.6 1 − 8 u , (5) dt 4 + Ma8 which can be thought of e.g. an optogenetic signal causing a direct induction of cyclin A expression. At high levels of cyclin A the effect of the optogenetic input rapidly declines mimicking saturation of the amount of cyclin A. Strength of this input is tuned by the factor which is essential in order to keep the oscillators in a neighborhood W of the limit cycle γ. For the calculation of the input u(t) the phase of each agent has to be determined. This process follows
is a local canonical model for such an oscillator. The phase θ(t) of an oscillator proceeds with constant angular velocity ω = 2π T and an additional term where the input u acts on the phase via the phase response curve Z(θ). The phase response curve (PRC) describes the magnitude of phase changes after perturbing an oscillatory system. Reduced phase models and the concept of PRCs are based on the notion of isochrons introduced by Winfree (1974) (Fig. 1). Therein, the set of all initial conditions x(0) ∈ Rn of which the solution x(t) approaches the solution ξ(t), with ξ(0) ∈ γ is called an isochron of ξ(0). The flow maps isochrons to isochrons and Guckenheimer (1975) showed, that there always exists a neighborhood W of a limit cycle that is invariantly foliated by isochrons. The function ψ2 : W → γ maps a point in the neighborhood x ∈ Mξ ⊂ W to the generator of its isochron ξ ∈ γ. Furthermore, the periodic orbit of an oscillator is homeomorphic to the unit circle and the function ψ1 : γ → S1 maps the solution ξ(t) with ξ(0) ∈ γ to the solution of Eq. (2). The function ψ : W → S1 is a composition of ψ1 and ψ2 , mapping x(t) ∈ W uniquely to its corresponding phase θ(t) of the reduced phase model (Fig. 1). Synchrony of a population of oscillators can be measured by means of the first circular moment which is defined as 2π N ˆ 1 (t) = m1 (t) = ei θk (t) , or m ei θ ρ(θ, t) d θ , (3) k=1
45
0
for individual agents θk (t) or the distribution of agents on the unit circle ρ(θ, t), respectively. In the complex number m1 = mˆ1 = r ei φ , φ ∈ S1 is the mean phase and 45
IFAC FOSBE 2018 46 Chicago, Illinois, USA, August 5-8, 2018
Karsten Kuritz et al. / IFAC PapersOnLine 51-19 (2018) 44–47
• Upon cell division, heterogeneity is introduced by sampling the initial condition of the two daughter cells from a logarithmic normal distribution with sampling variance σ and the mean equal to the state values at the limit cycle γ at cell division. • The total number of cells is kept constant by removing a randomly determined cell from the list cells whenever a cell divides. A snapshot measurement is taken at the end of each simulation loop with a simulation time of tˆend = 1 h from which the input for the next loop is calculated.
(a)
(c)
3.0
cyclin B/Cdk1
2.5
G2/M ↖ 2.0
↙ M/G1
1.5 1.0
S/G2 ↗
0.5 0.0 0
1
2
cyclin A/Cdk2
3
d Ma dt d Mb dt d Me dt d E2F dt
Ma Kda + Ma Mb =vsb · Ma − Vdb · Cdc20 · Kdb + Mb Me =vse · E2F − Vde · Ma · Kde + Me E2Ftot − E2F =V1e2 f · K1e2 f + (E2Ftot − E2F)
=vsa · E2F − Vda · Cdc20 ·
d Cdc20 Cdc20tot − Cdc20 =V1dc20 · Mb · dt K1cdc20 + (Cdc20tot − Cdc20) − V2cdc20 ·
Cdc20 K2cdc20 + Cdc20
Fig. 2. (a) Scheme of the Cdk network driving the mammalian cell cycle (G´erard and Goldbeter, 2011). (b) Limit cycle oscillations shown as a projection onto the cyclin A/Cdk2 versus cyclin B/Cdk1 phase plane. (c) Corresponding system of ODEs, with complexes cyclin A/Cdk2, cyclin B/Cdk1 and cyclin E/Cdk2 abbreviated by Ma, Mb and Me, respectively.
Fig. 3. Scheme of the first layer of the simulation framework (adapted from Imig et al. (2015)).
the theory of reduced phase models, where ψ2 : W → γ maps a point x(t) ∈ W to the generator of its isochron ξ(t) ∈ γ (Fig. 1). However, due to the high complexity of the notion of isochrons we used a simplified mapping where ψˆ 2 maps the point x(t) to its closest point on the ˆ limit cycle ξ(t) ∈ γ. The relation ψ1 : γ → S1 which maps a point on the limit cycle γ to its corresponding position on the unit circle S1 which is proportional to cell age can easily be obtained by inversion of the simulated trajectory on the limit cycle ξ−1 (γ) = θ/ω or by employing ergodic principles as described in Kuritz et al. (2017). The composition ψ = ψ1 ◦ ψˆ 2 of the two mappings now maps the point x(t) ∈ W uniquely to its phase θ(t) ∈ S1 on the unit circle. The PRC Z was obtained by solving the appropriate adjoint equation using the dynamic modeling program XPPAUT (Malkin, 1956; Ermentrout, 2002). The distribution of cells on the unit circle was calculated by kernel density estimation from phase values θk using a von Mises distribution as kernel with manually tuned bandwidth. The input is calculated from the PRC and the distribution of cells on the unit circle by Eq. (4). The simulation framework is based on the individual-based population model by Imig et al. (2015) and was slightly adapted, resulting in the following framework (Fig. 3): • The simulation runs from time-points t0 until tend . • A list cells contains all initial mother-cells at t0 . • The cell with the shortest simulated time t0 is simulated until a specified time tˆend , or until the cell complies with a division condition.
4. SIMULATION RESULTS The cell population was initialized by applying the inverse : S1 → γ on samples drawn from a von mapping ψ−1 1 Mises distribution with circular moment m1 (t0 ) = 0.3. Simulations with varying input factor and sampling variance σ were run for 500 hours. Exemplary results obtained for sampling variance of 0.13 and 0.18 and input strength varying between 0 and 0.09 illustrate the different effects of the controller on the cell population (Fig. 4).
length first moment
(a)
σ = 0.13
1 0.8 0.6 0.4 0.2 0
(b)
length first moment
(b)
0
100 200 300 400 500
time in h
σ = 0.18
1 0.8 0.6 0.4 0.2 0
0
100 200 300 400 500
time in h
Fig. 4. Absolute value of the first circular moment for sampling variances σ = 0.13 (a) and σ = 0.18 (b). For < 0.02, the input is too small and the norm of the first circular moment remains constant or decreases, indicating a conversion towards a steady-state distribution. Increasing the input strength to = 0.05 results in stable 46
IFAC FOSBE 2018 Chicago, Illinois, USA, August 5-8, 2018
Karsten Kuritz et al. / IFAC PapersOnLine 51-19 (2018) 44–47
that the proposed control strategy is robust such that synchronization might also be achieved under stochastic dynamics. Typically synchronization of a cell population is achieved by environmental or chemical treatments. This procedure is often accompanied by complex side effects and the population loses synchrony rapidly. The presented controller can be used to culture cell populations which are synchronized while hardly affecting individual cells in their natural cell cycle. REFERENCES
oscillations with increasing length of the first moment. However, input strengths ≥ 0.09 are too strong causing the cells to diverge from their original limit cycle which results in loss of oscillation of parts of the population. Cell cycle arrest of a fraction of the population results in a very large amplitude of the length of the first moment. When the oscillating subpopulation passes by the arrested cells |m1 | increases. When the cells are on opposite sites of the cell cycle |m1 | is small. Such a behaviour can be seen for = 0.09 at the early time points with σ = 0.13 and throughout the whole simulation with σ = 0.18 (Fig. 4). Input strength and sampling variance were varied over a sufficiently large range, as illustrated in Fig. 5, to cover relevant outcomes and define regions where: 1) the input predominates over noise and the population synchronizes, 2) the population diverges and noise introduced during cell division predominates, or 3) the input is too strong causing loss of oscillations. The input strength which is necessary to synchronize the population increases with the variance until it reaches a point where the input causes loss of oscillations. The border of this region expands at higher variances to lower input strengths. A possible explanation for this observation is that higher variances during resampling result in larger deviations from the limit cycle during cell division such that lower inputs suffice to interrupt stable oscillations. The decline rate of |m1 | in the simulation with zero input and sampling variances of σ = 0.2 is comparable to a cell cycle length distribution with mean 24 h and variance 2 h. Such a distribution is in good accordance with the experimentally observed cell cycle lengths of mammalian cell cultures.
Charvin, G., Cross, F.R., and Siggia, E.D. (2009). Forced periodic expression of G1 cyclins phase-locks the budding yeast cell cycle. Proceedings of the National Academy of Sciences of the United States of America, 106, 6632–6637. doi:10.1073/pnas.0809227106. Ermentrout, B. (2002). Simulating, Analyzing, and Animating Dynamical Systems. Society for Industrial and Applied Mathematics. G´erard, C. and Goldbeter, A. (2011). A skeleton model for the network of cyclin-dependent kinases driving the mammalian cell cycle. Interface Focus, 1(1), 24–35. Guckenheimer, J. (1975). Isochrons and phaseless sets. J. Math. Biol., 1(3), 259–273. Hanahan, D. and Weinberg, R.A. (2011). Hallmarks of cancer: the next generation. Cell, 144(5), 646–74. Hoppensteadt, F.C. and Izhikevich, E.M. (1997). Weakly Connected Neural Networks, volume 126 of Applied Mathematical Sciences. Springer New York. Huh, D. and Paulsson, J. (2011). Non-genetic heterogeneity from stochastic partitioning at cell division. Nature Genetics, 43(2), 95–100. doi:10.1038/ng.729. Imig, D., Pollak, N., Strecker, T., Scheurich, P., Allgower, ¨ F., and Waldherr, S. (2015). An individual-based simulation framework for dynamic, heterogeneous cell populations during extrinsic stimulations. J. Coupled Syst. Multiscale Dyn., 3(2), 143–155. Kuramoto, Y. (1984). Chemical Oscillations, Waves, and Turbulence, volume 19 of Springer Series in Synergetics. Springer Berlin Heidelberg. Kuritz, K., Stohr, D., Pollak, N., and Allgower, F. (2017). ¨ ¨ On the relationship between cell cycle analysis with ergodic principles and age-structured cell population models. J. Theor. Biol., 414(7), 91–102. Kuritz, K., Halter, W., and Allgower, F. (2018). Passivity¨ based ensemble control for cell cycle synchronization. In R. Tempo, S. Yurkovich, and P. Misra (eds.), Emerg. Appl. Control Syst. Theory. Springer International Publishing, 1 edition. Malkin, I.G. (1956). Some Problems in Nonlinear Oscillation Theory. Gostexizdat, Moscow. Milias-Argeitis, A., Rullan, M., Aoki, S.K., Buchmann, P., and Khammash, M. (2016). Automated optogenetic feedback control for precise and robust regulation of gene expression and cell growth. Nat. Commun., 7, 12546. Montenbruck, J.M., Burger, M., and Allgower, F. (2015). ¨ ¨ Practical synchronization with diffusive couplings. Automatica, 53, 235–243. Winfree, A.T. (1974). Patterns of phase compromise in biological cycles. J. Math. Biol., 1(1), 73–93. Zhivotovsky, B. and Orrenius, S. (2010). Cell cycle and cell death in disease: past, present and future. J. Intern. Med., 268(5), 395–409.
oscilla
division noise predominates
tions
0.2 0.15 population synchronizes 0.1
0
0.02
0.04 0.06 input strength
0.08
upted
interr
resampling variance σ
0.3 0.25
47
0.1
Fig. 5. Parameter regions where we observe no synchronization due to predominant division noise, successful synchronization or loss of oscillations. 5. CONCLUSION In this paper we describe a proof-of-concept for the recently derived ensemble control approach to synchronize a cell population within the cell cycle (Kuritz et al., 2018). We expand our previous work which was derived for an idealized behavior of the cell population to more realistic conditions. To this end we used individual-based simulation framework where heterogeneity is introduced by resampling of states upon cell division (Imig et al., 2015). The simulation results indicate the existence of a suitable range in which the input strength can be varied while still achieving synchrony without interrupting the oscillations. Heterogeneity in the population can also be introduced by stochastic dynamics. We hypothesize 47