I DENTIFYING DYNAMIC R EGULATION WITH A DVERSARIALS URROGATES

1

Ron Teichner1,2 , Naama Brenner3,2 , and Ron Meir1,2
Viterbi Dept. of Electrical & Computer Eng., Technion, Israel Institute of Technology
2
Network Biology Research Lab, Technion, Israel Institute of Technology
3
Dept. of Chemical Eng., Technion, Israel Institute of Technology

A BSTRACT
Homeostasis, the ability to maintain a stable internal environment in the face of perturbations, isessential for the functioning of living systems. Given observations of a system, or even a detailed
model of one, it is both valuable and extremely challenging to extract the control objectives of thehomeostatic mechanisms. Lacking a clear separation between plant and controller, frameworks suchas inverse optimal control and inverse reinforcement learning are unable to identify the homeostaticmechanisms. A recently developed data-driven algorithm, Identifying Regulation with Adversarial
Surrogates (IRAS), detects highly regulated or conserved quantities as the solution of a min-max
optimization scheme that automates classical surrogate data methods. Yet, the definition of homeostasis as regulation within narrow limits is too strict for biological systems which show sustained
oscillations such as circadian rhythms. In this work, we introduce Identifying Dynamic Regulation
with Adversarial Surrogates (IDRAS), a generalization of the IRAS algorithm, capable of identifying
control objectives that are regulated with respect to a dynamical reference value. We test the algorithm
on simulation data from realistic biological models and benchmark physical systems, demonstrating
excellent empirical results.
Keywords model discovery · control · adversarial · machine learning

1

Introduction

Living systems maintain stability against internal and external perturbations, a phenomenon known as homeostasis[7, 28, 31]. This is a ubiquitous central pillar across all scales of biological organization, such as molecular circuits,
physiological functions, and population dynamics. Failure of homeostatic control is associated with diseases including
diabetes, autoimmunity, and obesity [31]. It is therefore vital to identify the regulated variables that the system aims
to maintain at a stable setpoint. Unlike simple human-made systems, where known pre-selected variables are undercontrol, biological systems are characterized by multiple coupled control loops as well as other dynamic structures[7]. A fundamental difference is that they are not divided to separate ‘plant’ and ‘controller’ entities, as is commonly
assumed in control theory, but rather make up a complex network of interactions. In such a network regulated variables
may be maintained at a stable setpoint, as in the classic example of the regulation of blood glucose concentration, atightly regulated biological process in which the rates of glycolysis and gluconeogenesis are flexible variables [31].
A particularly interesting aspect is that, under certain conditions, biological homeostatic controllers may become
oscillatory [48]. Indeed, while the occurrence of oscillations is generally avoided in control engineering, oscillatory
behavior is ubiquitously found in natural systems [32, 17, 35]. A well documented example is the Baroreflex controlof the cardiovascular system where slow blood-pressure and heart-rate oscillations are observed [15]. The biologicalcomplexity makes it challenging to identify the regulated variables that the system actively maintains in the vicinity ofan attractor and, as we further explain, dedicated algorithmic tools are required.
Biological systems are commonly modeled by a set of dynamical equations where internal variables and control
variables are not clearly separated. Therefore, the controlled objective and the control signal in a biological systemmodel are implicit. Although control theory typically assumes the existence of a separate plant and controller (seeInternational Conference on Scientific Computing and Machine Learning 2024 (SCML2024)

SCML2024 - Identifying regulation algorithm

Chapter 1.2 in [6]), many theoretical results and analysis tools do not require such a separation [12]. Identifying
Regulation with Adversarial Surrogates (IRAS), a recently published algorithm directly addresses the plant-controller
coupling, and is thus dedicated to control identification in biological systems. IRAS identifies variables that are
regulated within narrow limits and is thus incompatible with biological systems which show sustained oscillations suchas circadian rhythms [38]. In this work we introduce Identifying Dynamic Regulation with Adversarial Surrogates(IDRAS), a fundamental extension of the IRAS algorithm, capable of identifying control objectives that are regulatedwith respect to a dynamic reference value. We test the algorithm on synthetic data of kinetic protein interactions,bacterial growth and division and on purely physical systems demonstrating excellent empirical results. See Appendix
Section C for a review of related work.
The IRAS algorithm assumes that it is in general a combination of the observed variables which is internally controlled.
Therefore, given a set of measurements over time, z(t) ∈ Rm , IRAS finds a combination g(·) of the measured variables,that is maintained around a setpoint,g(z(t)) ≈ cset ,(1)
where we can assume cset = 0 without loss of generality. In this work, we consider systems that are regulating some
unknown combination of the observables over an unknown dynamically changing reference value,g(z(t)) ≈ ce (t).

(2)

By adding the subscript e we emphasize that the behavior of ce (t) can be environment dependent. For example,
circadian rhythms are regulated on molecular and physiological levels, but control parameters are entrained to the
daylight cycle [39, 16]. Here we do not assume any specific form neither for the combination g() nor for the process
ce (). For illustration, a general model that might yield observations that satisfy g(z(t)) = ce (t) is,due = fue (ue )dt + dWue ,dx = fx (x, ue )dt + dWx ,zk = h(x(tk )) + ηzk ,(3)

where x ∈ Rn is the state of the observed system, ue ∈ Rnu is an input from the environment, zk ∈ Rm is theobservation sampled at time tk , dWue and dWx are Wiener processes and ηzk is a measurement noise.
In general, any system identification procedure that does not decouple g() from ce () will yield results that do not
generalize to other environments where ce′ (t) has a different behavior. It is therefore vital to separately identify theregulated combination from the dynamics of the reference. In the next Section we present IDRAS, a purely data-driven
algorithm that simultaneously learns a combination g() and the dynamic process that it follows, ce ().

2

Algorithm development

We are interested in identifying empirically, from a setof measurements, a variable combination which tightly
follows a dynamically changing reference value, where
both are unknown. The combination could represent aninternal quantity of high importance to the system, andthe reference value could reflect temporal trends in theenvironment. Given observations zk ∈ Rm at discrete
times k ∈ {1, 2, . . . , N }, we search for gθ : Rm → R, afunction parameterized by θ such that,ck = gθ (zk ),(4)

where the resulting time-series ck is the (learned) reference. We assume that this time-series follows somedynamics and it is therefore possible to learn a filter1 ,
namely a predictor of the value of ck at time k based on
its previous values, formally given byĉk = Fϕ,ω,δ (ck−T :k−1 ).

(5)

Figure 1: IDRAS algorithm outline. The observed time-series

z is permuted to create the unconstrained series z ∗ . The shuffle
player, only exposed to the 1D projections under EΩ (), sets theresampling function ζ() used to resample z̃ from z ∗ such thatthe distributions under the projection are identical, (8). Thenthe combination player, given z and z̃, updates the parametersΩ towards minimizing the CR, (7). The block EΩ , based on (6)replaces the function g() in Fig. 3 in [45].

The structure of the learned filter F , and the meaning of its parameters, is detailed at the end of the present Section, andck−T :k−1 = [ck−T , ck−T +1 , . . . , ck−1 ] with T a hyper-parameter. The filtering error is,∆

ek = ck − ĉk = EΩ (zk−T :k ),1

Formally (5) is a 1-step predictor, but, following [5], we refer to it as a filter.

2(6)

SCML2024 - Identifying regulation algorithm

and is a function of the parameters Ω = [θ, ϕ, ω, δ]. Our goal is to learn both gθ () and Fϕ,ω,δ () such that the error issmall, namely ek ≈ 0 for all k. We note that a straightforward optimization yields the trivial pair, g() ≡ 0, F () ≡ 0.
The problem formulated is in fact a generalization of problem (1). It is immediate to verify that by setting Fϕ,ω,δ (·) ≡ 0,
we obtain EΩ (zk−T :k ) = gθ (zk ), ek = ck and thus we optimize to find a combination gθ (·) such that ck ≈ 0, namely
regulated to a narrow range around a fixed value.
This simpler problem was addressed in [46] by utilizing comparisons between original time-series and random (thoughpossibly constrained) shuffled time-series. A regulated combination is presumably composed of components thatco-vary to compensate and buffer perturbations. Therefore, shuffling the temporal order of each component separately
would ruin these co-variations and greatly increase the variance of the combination over time. To quantify this notion, a“coefficient of regulation" (CR) is defined as the ratio between the standard deviations of the combination in the original
data ck and in the shuffled data c˜k . In [46] it was found that straightforward optimization of this measure is insufficientto escape trivial solutions and artefacts. Rather, a two-player algorithm was constructed which alternates between
minimizing the CR and modifying the shuffled time-series. This algorithm terminates when the CR can no longer beminimized, while at the same time reproducing general geometric features of the data in parameter space. Relying on
similar concepts, the next Section introduces the generalization of this algorithm to time-varying regulation set-points.
Identifying Dynamic Regulation with Adversarial Surrogates
IDRAS is an iterative algorithm consisting of two competing players, a generalization of the IRAS algorithm. The firstplayer aims to minimize the CR, a measure of invariance (see [46] Section 1), while the second player aims to render
the task of the first player more difficult by forcing it to extract information about the temporal structure of the data,
which is absent from time-shuffled “surrogate” data.
Combination player

At iteration i, the first player sets the parameters Ω to minimize the CR,Ω(i) = argmin σ(eT +1:N )/σ(ẽ(i)T +1:N ),

(7)

Ω
(i)

(i)

where σ(eT +1:N ) and σ(ẽT +1:N ) are the standard-deviations of the time-series ek and ẽk respectively. The time-series(i)
ek is the filtering error on the original data and the time-series ẽk is the filtering error on surrogate data, the outcome of(i−1)

a resampling procedure using a resampling function ζ() that was set by the second player in the previous iteration.

(i−1)
(i)
(i)
(i)
(i)
(i)
(i)
Formally, ek = EΩ (zk−T :k ), ẽk = EΩ (zk−T :k−1 z̃k ), and z̃k ∼ Pz̃ζ(z̃k | zk−T :k−1 ) = Pz (z̃k )ζ (i−1) (ẽk ),where Pz (z) is the probability of observing a measurement z at some random time.

Optimizing (7), the combination player searches for a combination gθ (·), whose filtering error (which is the 1-step(i)
prediction error) σ(eT +1:N ) is small w.r.t. the error in predicting a sample from a random point in time, σ(ẽT +1:N ).
This encourages the combination player to find a useful (g, F ) pair, representing a meaningful underlying quantity inthe observed system. We note that for ζ(·) ≡ 1 the time-series z̃k is a naïve random permutation of the time-series zk .
Optimizing (7) w.r.t. such unconstrained permutation leads to artifacts, see Appendix Section A for elaboration. To
avoid these artifacts the second player constrains the time-series z̃k by setting the resampling function ζ(·).
Shuffle player The second player uses the current solution of the combination player, Ω(i) , to create a new shuffled(i+1)
time-series z̃k, which resembles the statistical structure of the data under the 1D projection EΩ(i) . Formally, thiscorresponds to the selection of a resampling function ζ (i) (·) that minimizes the distributional distance,ζ (i) = argmin D(eT +1:N , ẽT +1:N ),ζ

(8)where ek = EΩ(i) (zk−T :k ), ẽk = EΩ(i) (zk−T :k−1 , z̃k ), and z̃k ∼ Pz̃ζ (z̃k | zk−T :k−1 ) = Pz (z̃k )ζ(ẽk ).
Lemma 1. (Shuffle player’s optimal solution.) The shuffle player, who only has access to the error time-series, can solve∗
optimization (8) and obtain D = 0 by choosing ζ (EΩ(i) (·)) = Pz(T ) (·)/Pz(T(·), where Pz(T ) (zk−T :k ) is the probability)∗ (zk−T :k−1 , z) = Pzof observing the sequence zk−T :k (at a random time k) and Pz(T(zk−T :k−1 )Pz (z) is the(T −1))
probability of observing the sequence zk−T :k−1 followed by a time-random observation, sampled from the naïvly
permuted time-series.
The two players inform each other of their current step results, and the process continues iteratively until the combination
player can no longer decrease the CR in (7). We refer to this algorithm as IDRAS and depict its outline in Figure 1. The
filter block (5), the proof of Lemma 1 and a performance assessment metric are given in the Appendix Section A.
3

SCML2024 - Identifying regulation algorithm

3

Validation

After presenting the construction of IDRAS, we seek
to validate it on datasets with a known control objective, so that the quality of the results can be assessed.
We chose two validation examples of biological models: a kinetic model of protein interactions and a modelof bacterial life cycle. Additionally, to demonstrate theefficiency of our algorithm in studying physical systems, we validate the algorithm on a dataset that serves
for benchmarking machine-learning algorithms. We
briefly depict here the first examples and leave a detailed explanation of all examples to the Appendix Section B. Code reproducing all examples is available athttps://github.com/RonTeichner/IRAS.
A kinetic model of interactions
A kinetic model describes regulatory interactions in the
production of two proteins incorporating a feedback loop. Figure 2: (a) Illustration of the closed-loop model (9). The
In the considered model (inspired by [22]), the total mRNA molecule M induces the production of proteins S andamount of two proteins P and S, namely P + S, is P and receives a negative feedback of their sum. (b) The IRAS
controlled by M , the mRNA molecule that is produced algorithm does not converge, due to the lack of a regulated constantthrough a process of constitutive transcription at a rate combination. The combination P +S (black) along with the output
K and degraded with first order kinetics at a rate γM . of IRAS (dashed red). (c) IDRAS algorithm captures the control
objective and its oscillating trend.
The model is under constant perturbations to the protein
expression rate K. The model is described by the differential equations
1
K(t) = K0 (1 + 0.5 cos(2π t + ϕK )),dM = (K(t) − f (P + S) − γM M )dt,
τK
(9)
dP = (kP M − γP P )dt + ηP dWP ,dS = (kS M − γS S)dt + ηS dWS ,where the mRNA M and the two proteins P, S are linked in a feedback loop. Small changes in S or P , modeled by the
an increments of the Wiener processes dWP and dWS , induce swift and sharp changes in the transcription of M andmaintain P +S around a reference level−1
γM γP γS∆c∗ (t) = (P + S)(t) = f +
K(t),
(10)
kP γS + kS γP
which is reflected in a high negative correlation between S and P . In (9) we identify the components of the general model(3): ue = [K, M ], x = [P, S]. Our observations contain the levels of the two proteins P and S (an observation taken inlabs by measuring fluorescence intensity of constitutively expressed proteins [11]), zk = [P (t = k/fs ), S(t = k/fs )],where k ∈ [1, . . . , N ] and fs is the sampling rate. Figure 4a illustrates the kinetic interactions model, parameter values
are listed in the Appendix.
We ran both IRAS and IDRAS algorithms in search of the control objective. Figure 4b depicts the output (dashed red)of gθ (P, S) trained by the IRAS algorithm. Due to the lack of a regulated constant combination, IRAS did not converge.
When running IDRAS a combination that tightly follows a dynamic reference was found despite its oscillating nature,
as depicted in Figure 4c. See the Appendix for performance assessment metric scores.

4

Discussion

Biological systems are often “reversed engineered" to build mathematical models based on observed data [34, 41, 13,8, 42, 25, 10]. Detecting regulatory processes is a technically challenging problem with many potential applications[1, 36, 2, 18, 20, 19, 21]. Recently an algorithm named IRAS was introduced [46] that receives as input raw dynamic
measurements and provides combinations of the observables that are maximally conserved across time. Here we
presented IDRAS, a generalization of IRAS, capable of identifying control objectives that are regulated w.r.t. adynamical reference value. This algorithm allows the identification of control objectives in biological systems whichshow sustained oscillations, such as circadian rhythms. We provide validation in three distinct realistic examples
demonstrating excellent empirical results.
4

SCML2024 - Identifying regulation algorithm

References[1] Abdullahi Rashid Adem and Chaudry Masood Khalique. Symmetry reductions, exact solutions and conservation
laws of a new coupled kdv system. Communications in Nonlinear Science and Numerical Simulation, 17(9):3465–
3475, 2012.
[2] Abdullahi Rashid Adem and Ben Muatjetjeja. Conservation laws and exact solutions for a 2d zakharov–kuznetsov
equation. Applied Mathematics Letters, 48:109–117, 2015.
[3] Ferran Alet, Dylan Doblar, Allan Zhou, Josh Tenenbaum, Kenji Kawaguchi, and Chelsea Finn. Noether networks:
meta-learning useful conserved quantities. Advances in Neural Information Processing Systems, 34:16384–16397,
2021.
[4] Ariel Amir. Cell size regulation in bacteria. Physical review letters, 112(20):208102, 2014.
[5] Brian DO Anderson and John B Moore. Optimal filtering. Courier Corporation, 2012.
[6] Karl Johan Åström and Richard M Murray. Feedback systems: an introduction for scientists and engineers.
Princeton university press, 2021.
[7] George E Billman. Homeostasis: the underappreciated and far too often ignored central organizing principle ofphysiology. Frontiers in Physiology, page 200, 2020.
[8] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse
identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932–
3937, 2016.
[9] Michael B Chang, Tomer Ullman, Antonio Torralba, and Joshua B Tenenbaum. A compositional object-based
approach to learning physical dynamics. arXiv preprint arXiv:1612.00341, 2016.
[10] Boyuan Chen, Kuang Huang, Sunand Raghupathi, Ishaan Chandratreya, Qiang Du, and Hod Lipson. Automateddiscovery of fundamental variables hidden in experimental data. Nature Computational Science, 2(7):433–442,
2022.
[11] Yong-Joon Choi and Kazuaki Sawada. Physical sensors: Fluorescence sensors. In Roger Narayan, editor,
Encyclopedia of Sensors and Biosensors (First Edition), pages 1–19. Elsevier, Oxford, first edition edition, 2023.
[12] Carlo Cosentino and Declan Bates. Feedback control in systems biology. Crc Press, 2011.
[13] Bryan C Daniels and Ilya Nemenman. Automated adaptive inference of phenomenological dynamical models.
Nature communications, 6(1):1–8, 2015.
[14] Filipe de Avila Belbute-Peres, Kevin Smith, Kelsey Allen, Josh Tenenbaum, and J Zico Kolter. End-to-end
differentiable physics for learning and control. Advances in neural information processing systems, 31, 2018.
[15] M Di Rienzo, G Parati, A Radaelli, and P Castiglioni. Baroreflex contribution to blood pressure and heart rateoscillations: time scales, time-variant characteristics and nonlinearities. Philosophical Transactions of the Royal
Society A: Mathematical, Physical and Engineering Sciences, 367(1892):1301–1318, 2009.
[16] Jeanne F Duffy and Kenneth P Wright Jr. Entrainment of the human circadian system by light. Journal of
biological rhythms, 20(4):326–338, 2005.
[17] Jay C Dunlap, Jennifer J Loros, and Patricia J DeCoursey. Chronobiology: biological timekeeping. Sinauer
Associates, 2004.
[18] OH El-Kalaawy. Variational principle, conservation laws and exact solutions for dust ion acoustic shock waves
modeling modified burger equation. Computers & Mathematics with Applications, 72(4):1031–1041, 2016.
[19] OH El-Kalaawy. New: Variational principle–exact solutions and conservation laws for modified ion-acoustic
shock waves and double layers with electron degenerate in plasma. Physics of Plasmas, 24(3):032308, 2017.
[20] OH El-Kalaawy. Modulational instability: conservation laws and bright soliton solution of ion-acoustic waves inelectron-positron-ion-dust plasmas. The European Physical Journal Plus, 133(2):1–12, 2018.
[21] OH El-Kalaawy, SM Moawad, and Shrouk Wael. Stability: Conservation laws, painlevé analysis and exact
solutions for s-kp equation in coupled dusty plasma. Results in physics, 7:934–946, 2017.
[22] Hana El-Samad. Biological feedback control—respect the loops. Cell Systems, 12(6):477–487, 2021.
[23] Richard P Feynman, Robert B Leighton, and Matthew Sands. The feynman lectures on physics; vol. i. American
Journal of Physics, 33(9):750–752, 1965.
[24] Jacopo Grilli, Matteo Osella, Andrew S Kennard, and Marco Cosentino Lagomarsino. Relevant parameters inmodels of cell division control. Physical Review E, 95(3):032411, 2017.
5

SCML2024 - Identifying regulation algorithm[25] Adam Haber and Elad Schneidman. Learning the architectural features that predict functional similarity of neuralnetworks. Physical Review X, 12(2):021051, 2022.
[26] Jessica B Hamrick, Kelsey R Allen, Victor Bapst, Tina Zhu, Kevin R McKee, Joshua B Tenenbaum, andPeter W Battaglia. Relational inductive bias for physical construction in humans and machines. arXiv preprint
arXiv:1806.01203, 2018.
[27] Po-Yi Ho, Jie Lin, and Ariel Amir. Modeling cell size regulation: From single-cell-level statistics to molecularmechanisms and population-level effects. Annual review of biophysics, 47:251–271, 2018.
[28] Victoria Hsiao, Anandh Swaminathan, and Richard M Murray. Control theory for synthetic biology: recent
advances in system characterization, control design, and controller implementation for synthetic biology. IEEE
Control Systems Magazine, 38(3):32–62, 2018.
[29] Suckjoon Jun and Sattar Taheri-Araghi. Cell-size maintenance: universal strategy revealed. Trends in microbiology,
23(1):4–6, 2015.
[30] Patrick Kidger, James Foster, Xuechen Li, Harald Oberhauser, and Terry Lyons. Neural SDEs as InfiniteDimensional GANs. International Conference on Machine Learning, 2021.
[31] Maya E Kotas and Ruslan Medzhitov. Homeostasis, inflammation, and disease susceptibility. Cell, 160(5):816–827,
2015.
[32] Dorothy T Krieger. The clocks that time us:(physiology of the circadian timing system). Psychosomatic Medicine,44(6):559–560, 1982.
[33] Xuechen Li, Ting-Kam Leonard Wong, Ricky T. Q. Chen, and David Duvenaud. Scalable gradients for stochastic
differential equations. International Conference on Artificial Intelligence and Statistics, 2020.
[34] Lennart Ljung. System identification. In Signal analysis and prediction, pages 163–173. Springer, 1998.
[35] David Lloyd and Ernest L Rossi. Ultradian rhythms in life processes: An inquiry into fundamental principles ofchronobiology and psychobiology. Springer Science & Business Media, 2012.
[36] Stanislav Yu Lukashchuk. Conservation laws for time-fractional subdiffusion and diffusion-wave equations.
Nonlinear Dynamics, 80(1):791–802, 2015.
[37] Liang Luo, Yang Bai, and Xiongfei Fu. Stochastic threshold in cell size control. Physical Review Research,
5(1):013173, 2023.
[38] Martin C Moore-Ede. Physiology of the circadian timing system: predictive versus reactive homeostasis. American
Journal of Physiology-Regulatory, Integrative and Comparative Physiology, 250(5):R737–R752, 1986.
[39] Colin S Pittendrigh and Dorothea H Minis. The entrainment of circadian oscillations by light and their role as
photoperiodic clocks. The American Naturalist, 98(902):261–294, 1964.
[40] Adam Santoro, David Raposo, David G Barrett, Mateusz Malinowski, Razvan Pascanu, Peter Battaglia, andTimothy Lillicrap. A simple neural network module for relational reasoning. Advances in neural information
processing systems, 30, 2017.
[41] Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324(5923):81–
85, 2009.
[42] Jingxiang Shen, Feng Liu, Yuhai Tu, and Chao Tang. Finding gene network topologies for given biological
function with recurrent neural network. Nature communications, 12(1):1–10, 2021.
[43] Kate Sidaway-Lee, Maria J Costa, David A Rand, Bärbel Finkenstadt, and Steven Penfield. Direct measurement oftranscription rates reveals multiple mechanisms for configuration of the arabidopsis ambient temperature response.
Genome biology, 15:1–10, 2014.
[44] Lee Susman, Maryam Kohram, Harsh Vashistha, Jeffrey T Nechleba, Hanna Salman, and Naama Brenner.
Individuality and slow dynamics in bacterial growth homeostasis. Proceedings of the National Academy of
Sciences, 115(25):E5679–E5687, 2018.
[45] Ron Teichner, Ron Meir, and Danny Eitan. Enhancing causal estimation through unlabeled offline data. arXiv
preprint arXiv:2202.07895, 2022.
[46] Ron Teichner, Aseel Shomar, Omri Barak, Naama Brenner, Shimon Marom, Ron Meir, and Danny Eytan. Identifying regulation with adversarial surrogates. Proceedings of the National Academy of Sciences,120(12):e2216805120, 2023.
[47] Joshua B Tenenbaum, Vin de Silva, and John C Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290(5500):2319–2323, 2000.
6

SCML2024 - Identifying regulation algorithm[48] Kristian Thorsen, Oleg Agafonov, Christina H Selstø, Ingunn W Jolma, Xiao Y Ni, Tormod Drengstig, and Peter
Ruoff. Robust concentration and frequency control in oscillatory homeostats. PLoS One, 9(9):e107766, 2014.
[49] Silviu-Marian Udrescu and Max Tegmark. Ai feynman: A physics-inspired method for symbolic regression.
Science Advances, 6(16):eaay2631, 2020.
[50] Nicholas Watters, Daniel Zoran, Theophane Weber, Peter Battaglia, Razvan Pascanu, and Andrea Tacchetti. Visualinteraction networks: Learning a physics simulator from video. Advances in neural information processing
systems, 30, 2017.
[51] Hong Wei Yang, Min Guo, and Hailun He. Conservation laws of space-time fractional mzk equation for rossby
solitary waves with complete coriolis force. International Journal of Nonlinear Sciences and Numerical Simulation,20(1):17–32, 2019.
[52] Aston Zhang, Zachary C. Lipton, Mu Li, and Alexander J. Smola. Dive into Deep Learning. 2020. https://d2l.ai.

A

Appendix A: Algorithm development

Disabling the shuffle player in the algorithm, that is, setting ζ(·) ≡ 1, the time-series z̃k is a naïve random permutationof the time-series zk . Optimizing (7) w.r.t such unconstrained permutation lead to artifacts. To illustrate consider that
one of the observables is merely a Wiener process. The variance of a 1-step prediction error is proportional to dt
while the variance of predicting a random sample increases with N , the length of the time-series. Therefore, without
constraining the permutation, the CR of a trivial combination that outputs the Wiener process approaches zero as N
increases. This leads to the identification of the Wiener process as the control objective, clearly an artifact. We refer thereader to Section 1.B in [46] for a proof and a detailed explanation.
Proof of Lemma 1. In what follows we prove that the resampling function
ζ (EΩ(i) (·)) =Pz(T ) (·)
,
∗ (·)
Pz(T
)

(11)

z̃k ∼ Pz̃ζ (z̃k | zk−T :k−1 ) = Pz (z̃k )ζ(ẽk ).

(12)

defined in 1 yields D = 0 in (8). Recall that,ek = EΩ(i) (zk−T :k ), ẽk = EΩ(i) (zk−T :k−1 , z̃k ),Based on (12) definePz̃(T ) (zk−T :k−1 , z̃k ) = Pz(T −1) (zk−T :k−1 )Pz̃ζ (z̃k | zk−T :k−1 ),the probability that the shuffle player will concatenate z̃k to the observed sequence zk−T :k−1 . In (8)-(12), the time-series
ek inherits its distribution via the 1D projection from the distribution Pz(T ) and the time-series ẽk inherits its distribution
from Pz̃(T ) . We will show that for the choice of the resampling function in (11), the two 1D projections coincide. Notethat
Pz̃(T ) (zk−T :k−1 , z̃k ) = Pz(T −1) (zk−T :k−1 )Pz̃ζ (z̃k | zk−T :k−1 )(i)= Pz(T −1) (zk−T :k−1 )Pz (z̃k )ζ(ẽk )(ii)∗ (zk−T :k−1 , z̃k )ζ(ẽk ),= Pz(T)(13)

where (i) is by the definition of Pz̃ζ in (12) and (ii) by the definition of Pz∗ (T ) in Lemma 1. The derivation in (13)implies that
∗ (ẽk )ζ(ẽk ),Pz̃(T ) (ẽk = EΩ(i) (zk−T :k−1 , z̃k )) = Pz(T)
Pz(ẽk ))
from which by substituting (11), ζ(ẽk ) = P (T∗(ẽk ) , it directly follows that Pz̃(T ) (ẽk ) = Pz(T ) (ẽk )
z

A.1

□.

(T )

Filter

The architecture of the filter within the block EΩ () in Fig. 1 has many degrees-of-freedom and can be chosen by theuser according to prior knowledge regarding the nature of the dynamic reference. To impose few constrains on the filter
F (·) it can be implemented by a fully connected deep neural-network. Dealing with biological systems, we assume
7

SCML2024 - Identifying regulation algorithm

Figure 3: Architecture of the EΩ (·) block (6) from Fig. 1. The filtering error ek is the difference between the learned+
reference value ck and its 1-step prediction ĉk . The filter Fϕ,ω,δ (·), (14), infers a latent-state yk−1
and time-advances it−to the state yk from which the estimation ĉk is decoded.
that the dynamics of the reference can be modeled by a continuous-time, time-invariant latent model (see [33, 30]
for details on integrating differential equations using neural-networks). The filter-block contains three parts: (i) An+
encoder eϕ : RT → Rny , that given T consecutive values of the reference, ck−T :k−1 , infers a latent state yk−1
∈ Rn yny
ny(with ny a user-defined hyper-parameter), (ii) A drift function wω : R → R describing the deterministic term+in the dynamics of the latent state that serves to time-advance the latent state yk−1
to yk− , (iii) An emission function−nydδ : R → R that decodes the 1-step predicted value ĉk from the latent state yk . The following set of equations
describe the filter-block (5) marked EΩ () in 1,
Z tk
+
+
yk−1
= eϕ (ck−T :k−1 ), yk− = yk−1+wω (y(t))dt, ĉk = dδ (yk− ).
(14)
tk−1

The block is depicted in Figure 3.
A.2

Performance assessment metric

IDRAS, an unsupervised learning algorithm, yields a time-series of a regulated quantity, ck , and the correspondingtime-series of 1-step predictions, ĉk . To assess whether a quantity that follows a dynamic reference value was indeed
found, we calculate the normalized prediction error energy ρ(c, ĉ).
Definition A.1. The normalized error energy for two time-series x1:N and x̂1:N is defined,PN¯k )2(x̄k − x̂
xk − E[x]ρ(x1:N , x̂1:N ) = k=1,, x̄k =
PN
2
σ(x)k=1 (x̄k )PN
PN
where E[x] = N1 k=1 xk and σ 2 (x) = N1 k=1 (xk − E[x])2 , and similary for x̂.
The normalized error energy ρ(c, ĉ) allows for an immediate assessment of the significance of the identified combination.
A low ρ value (ρ < 0.5) indicates that a combination that follows a dynamic reference was found. We remind thereader that the shuffle player in IDRAS guarantees that the algorithm won’t converge to trivial combinations such as
ck ≡ const.

B

Appendix B: Validation

In all examples, to assess the performance, we calculate the normalized error energy between the known control
objective c∗k and the output of IDRAS ck and consider ρ(c, c∗ ) < 0.1 as an excellent empirical result. We start byintroducing the architecture of the networks used throughout the examples. Code reproducing all examples is available
at https://github.com/RonTeichner/IRAS.
B.1

Architecture

All the parameterized functions used in Section 3 are feed-forward artificial neural networks with an input layer of
32 neurons, a hidden layer of 16 neurons and an output layer of appropriate dimension. The activation function of
8

SCML2024 - Identifying regulation algorithm

all neurons is Leaky-ReLU except for the output neuron whose activation function is the Sigmoid function. The
hyper-parameters are T = 10 and ny = 2. We train the networks as described in Algorithm 1 for a pre-defined fixed
number of epochs (150) using the Stochastic-Gradient-descent optimizer with a momentum value of 0.9 and a learning
rate of 0.01. See [52] for a detailed explanation of feed-forward neural networks, activation functions and optimizers.
B.2

A kinetic model of interactions

We first validate IDRAS on simulated data generated from a kinetic model that describes regulatory interactions in
the production of two proteins incorporating a feedback loop. In the considered model (inspired by [22]), the total
amount of two proteins P and S, namely P + S, is controlled by M , the mRNA molecule that is produced through aprocess of constitutive transcription at a rate K and degraded with first order kinetics at a rate γM . The model is underconstant perturbations to the protein expression rate K. These might be caused by ambient temperature dynamics thatwere found to have large passive effects on both mRNA synthesis and decay rates [43]. The model is described by thedifferential equations
1
t + ϕK )),τK
dM = (K(t) − f (P + S) − γM M )dt,
dP = (kP M − γP P )dt + ηP dWP ,dS = (kS M − γS S)dt + ηS dWS ,K(t) = K0 (1 + 0.5 cos(2π(15)

where the mRNA M and the two proteins P, S are linked in a feedback loop. Both P and S are positively affectedby M , with their steady-state values proportional to it. The concentration M in turn, is negatively affected by thesum P + S, with the strength of this negative feedback given by the rate constant f . We note that (15) is a model ofdimension 5, comprising 3 state variables and an oscillating input K that results from an underlying state-space modelof dimension at least 2 (a 1D model cannot produce oscillations).
Small changes in S or P , modeled by the an increments of the Wiener processes dWP and dWS , induce swift andsharp changes in the transcription of M and maintain P +S around a reference level∆c∗ (t) = (P + S)(t) =1K(t),γM γP γSf+
kP γS + kS γP(16)which is reflected in a high negative correlation between S and P [43]. In (15) we identify the components ofthe general model (3): ue = [K, M ], x = [P, S]. Our observations contain the levels of the two proteins P
and S (an observation taken in labs by measuring fluorescence intensity of constitutively expressed proteins [11]),zk = [P (t = k/fs ), S(t = k/fs )], where k ∈ [1, . . . , N ] and fs is the sampling rate. Figure 4a illustrates the kinetic
interactions model and SI Appendix section 1.2 lists the parameter values.
We ran both IRAS and IDRAS algorithms in search of the control objective. Figure 4b depicts the output (dashed red)of gθ (P, S) trained by the IRAS algorithm. Due to the lack of a regulated constant combination, IRAS did not converge
and scored ρIRAS (c, c∗ ) = 1.82 (see Definition A.1).
When running IDRAS, we have a simple measure to assess whether a quantity that follows a dynamic reference value wasfound, the measure of the normalized prediction error energy defined in Section A.2. A score of ρIDRAS (c, ĉ) = 0.148
was obtained, indicating that a combination that tightly follows a dynamic reference was found. The validation score isρIDRAS (c, c∗ ) = 0.012, indicating that the combination P + S was precisely found, despite its oscillating nature, as
depicted in Figure 4c.
Following are the exact details of the parameters used.
The timescale of the feedback loop τ(P +S) = 1/f = 0.0005(F = 2000) is much shorter than the timescale ofenvironmental influenced oscillations in K, τK = 0.2. We simulated a dataset of 100 observed systems sampled at rate
fs = 1 with Ks = Kp = 150, γS = γP = 70, γM = 80, f = 2000, K0 = 300 and dWP , dWS ∼ N (0, (0.5dt)2 ).
For each system ϕK ∼ U[0, π] and M (0), P (0), S(0) ∼ U[0.02, 0.1].
B.3

Bacterial life cycle

The next example we consider is a realistic biological model of bacterial life cycle. We shall focus here on bacterial
growth homeostasis, where growth and division proceed for many generations with significant variability and statistical
9

SCML2024 - Identifying regulation algorithm

Figure 4: (a) Illustration of the closed-loop model (15). The mRNA molecule M induces the production of proteins S
and P and receives a negative feedback of their sum. (b) The IRAS algorithm does not converge, due to the lack ofa regulated constant combination. The combination P + S (black) along with the output of IRAS (dashed red). (c)IDRAS algorithm captures the control objective and its oscillating trend.
stability. This is a problem with a long history and on which a large body of data is available. We apply our algorithm
to simulation data, which mimic experimental measurements but where the regulation is known, and show that IDRAS
detects the correct mode of regulation.
Most bacteria grow smoothly and divide abruptly, consistent with a threshold crossing by some division indicator at thesingle-cycle timescale [4, 24, 44]; the threshold itself can have dynamics over multiple cycles. Three types of divisionindicators corresponding to different regulation modes have been proposed: cell size (“sizer” control mechanism),
added size (“adder” mechanism) and elapsed time (“timer”) [29, 27]. In a recent paper, Luo et al. [37] demonstratethat commonly used heuristic tools uncover the correct mode of regulation only under very restricted conditions.
We reproduce these results and show that in contrast to the commonly used tools, IDRAS accurately identifies thecell-division mechanism.
The inter-generation threshold u(t) is modeled by a stochastic Ornstein-Uhlenbeck process,s
σ2µu − udt + 2 u dW
du =
τu
τu

(17)

where the term dW is an increment of a Wiener process [37]. On the single-cycle scale, the k th cell grows exponentially,x(t) = xk,b eαk (t−tk,b ) ,xk,b = x(tk−1,d )fk ,tk,b ≤ t ≤ tk,dtk,b = tk−1,dtk,d = argmint {t ≥ tk,b | xk,b e∆α(t−tk,b )(18)= u(t)},∆

where xk,b = x(tk,b ) is the birth size, xk,d = x(tk,d ) is the size at division, αk is the exponential growth rate and fk isthe division fraction. Here we simulate a sizer mechanism such that the cell divides when its size crosses the thresholdu(t). Figure 5a depicts a simulated lineage over time - the cell-size, x(t) (dashed-black), and the stochastic thresholdu(t) (blue). SI Appendix section 1.3 lists the parameter values.
In (17) and (18) we identify the components of the general model, (3), with ue = u and x = x. Our observations,derived from x(t), contain the initial size, growth rate and the cycle duration, zk = [xk,b , αk , Tk = tk,d − tk,b ] where
αk = T1k log( xdxb ). We note that our choice of the feature vector zk renders IDRAS’s task harder as now to correctly
detect the sizer mechanism the network has to learn the combination g() = xb eαT and not just g() = xd , in case xd
was an entry of zk .
We ran both IRAS and IDRAS algorithms in search of the division mechanism. We chose a realistic parameter set for
which heuristic identification methods based on data correlations fail to detect the control mechanism correctly [37].
Figure 5b depicts the output (dashed red) of gθ (xb , α, T ) trained by the IRAS algorithm vs the ground-truth output, thesizer division mechanism xb eαT . Here IRAS, optimizing objective (1), outputs a combination that results from fusing
together the division mechanism and the dynamic threshold u(t) (with ρIRAS (c, c∗ ) = 0.66), as we now explain. In ourmodel, the timescale of the threshold is significantly slower than the time scale of single-cells, τu ≫ T , such that the10

SCML2024 - Identifying regulation algorithm

Figure 5: (a) Bacteria life-cycles. The blue curve is the stochastic threshold process u(t), (17), and the black curveis the cell-size over time, (18), such that the cell divides once it reached the dynamic threshold.(b) IRAS does notdecouple the cell-size mechanism (black-curve) from the threshold trend c(t) and yields a combination (dashed-red)
that represents their mixture. (c) IDRAS captures the division-mechanism.
threshold is approximately constant along a single-cell cycle, u(t) ≈ uk for tk,b ≤ t ≤ tk,d . The birth size, about half
the size at division of the previous cell is thus xk,b ≈ 0.5uk , and the size at division is (by definition) xk,d ≈ uk . The
single-cell cycle regulated combination is therefore,
xk,b eαk Tk − uk ≈ xk,b eαk Tk − 2xk,b ≈ 0,a mixture of two terms. The first, xk,b eof the dynamic threshold process.

αk Tk

(19)

, represents the sizer mechanism while the second, −2xk,b , is an influence

IDRAS, optimizing objective (2), decouples the two by separately learning the threshold dynamics and a combination
which is regulated w.r.t. this dynamic reference value. Figure 5c depicts the precise identification of the sizer mechanismwith a validation score ρIDRAS (c, c∗ ) = 0.046. The normalized prediction error energy is ρIDRAS (c, ĉ) = 0.23, fitting the
expected value derived from the noise term in (17) up to 8% (2T̄/τc , where T̄ is the mean cell-cycle, about 25 minutes).
Following are the exact details of the parameters used.
The simulated dataset contains 100 lineages each consisting of 100 generations. In (18), µu = 1 [µm], σu = 0.1 [µm]and τu = 200 [min]; Initial conditions are u(0) ∼ N (µu , σu2 ) for each lineage. In Eq. 18, the growth rate is sampledfrom a Gamma-distribution, αn ∼ Γ(25, 9.4 · 10−4 ) and the division fraction is distributed fn ∼ N (0.5, 0.052 ). The
initial conditions for the first cell in each lineage are x0,b ∼ N (0.5, 0.052 ), t0,b = 0.
B.4

Identifying Complex Physical Equations

To further challenge the IDRAS algorithm, we evaluate it on a broad range of physics problems, taken from the seminal
Feynman Lectures on Physics [23], also used in the recently published “Feynman Symbolic Regression Database”(FSReD) [49]. We simulate physical equations incorporating oscillating terms and expect IDRAS to decouple theoscillating terms from the quantities they follow.
One such equation is equation II.6.15b, taken from Volume 2 in [23], which describes the electric field induced by adipole. The dipole, with a dipole moment Pd (the product of the charges and their separation) induces, at distance r andangle ωt, an electric field whose transverse component is Ef ,
3 Pd
Ef =
cos(ωt) sin(ωt).
4πϵ r3We simulate 100 time-series of length 100 in which the dipole spins with angular velocity ω. For each sample, theobservables (Ef , ϵ, Pd and r) are sampled uniformly in [1, 5] (although not physically valid, it is the way the FSReD
dataset is synthesized [49]). We expect IDRAS to identify the combination
4πϵr3gθ∗ (Ef , ϵ, Pd , r) = Ef
= cos(ωt) sin(ωt) = c∗ (t).
3Pd
To quantify the performance of IDRAS, we compute ρ(c, c∗ ) for each of three chosen examples from the fields ofelectric-field forces (I.12.11), non-linear responses (I.50.26), and electric dipoles II.6.15b. As shown in table 1, theagreement is excellent.
11

SCML2024 - Identifying regulation algorithm

Table 1: IDRAS captures physical relations

I.12.11
I.50.26
II.6.15b

C

Equationρ(c, c∗ )F = q ∗ (Ef + Bv sin(ωt))
x = x1 (cos(ωt) + α cos2 (ωt))
3 Pd
Ef = 4πϵr 3 cos(ωt) sin(ωt

0.007
0.022
0.059

Appendix C: Related work

Some dynamical systems obey conservation laws, where a combination of variables is constant along time trajectoriesof the system [51], as in (1). Identifying conservation laws from observed data is an active research field, [50, 40, 26,
14, 9, 47], and it is important to note that, even given the differential equations of the system, identifying a conservation
law analytically - or even proving its existence - is a difficult mathematical problem [1, 36, 2, 18, 20, 19, 21]. Methods
for deriving the differential equations of a system [8] leave the question of identifying the conservation law unsolved.
This problem is a sub-space of the one considered here, since the concept of conservation laws usually refers towell-defined dynamical systems with fixed parameter. Identifying biological regulation includes also the more generalcase where parameters are varying in time but still some combination is regulated around a set-point in the face of these
perturbations. Optimizing for conservation alone can lead to trivial quantities, such as predicting a constant g(z) = c
independent of z. In a recent paper, [3] refer to a non-trivial g(·) by the term useful conservation law. To obtain anon-trivial solution, [46] define a measure of invariance, the “Coefficient of Regulation" (CR) and an optimization
algorithm that obtains meaningful invariants in (1). Inferring meaningful dynamic regulatory processes, (2), renders thetask more difficult. Here we generalize IRAS to allow inferring meaningful dynamic regulatory process.

12