Time delays play a very important role in genetic regulatory systems. Gene regulation and signal transduction as a whole involves the synthesis and maturation of complex proteins. Their synthesis and transport takes a considerable amount of time, which introduces delays in the overall regulation chain. At a process level, metabolic time delays can be observed macroscopically by recognizing a certain time delay between substrate uptake and the corresponding biomass growth or product formation as in cultivations of *Saccharomyces cerevisiae* [1] or *Pichia pastoris* [2]. The nature of time delays in regulatory networks is twofold [3]. They are either related to a process that takes an intrinsic time to be accomplished, i.e. some reactions, such as translational or transcriptional reactions, take a significant amount of time to be completed, or as a consequence of the modelling approach used, i.e. lumping a sequence of events might lead to an apparent time delay.

The bottom-up systems biology approach for building dynamic network models can be too cumbersome due to their complex nature and lack of fundamental knowledge [4–7]. Typical limitations are the involvement of large scale kinetic models with poorly defined kinetic parameters, limited generalization capacity and their cost expansive development. In this paper we propose the use of mathematical hybrid semi-parametric systems as a cost effective alternative to model biochemical networks with intrinsic time delays, since it is not likely to know in advance which fundamental mechanisms cause such delays. Hybrid semi-parametric systems combine fundamental (parametric) biological constraints with more empirical data-based (nonparametric) constraints. Mechanistic and nonparametric models can therein be arranged in two possible ways: parallel or serial [4, 8–12]. In the serial structure, which has been the one applied in this study, the biological system dynamics are described by time differentials of classifying variables, while the unknown metabolic functions with intrinsic delays are handled by a nonparametric structure.

A general mathematical representation of delayed dynamics is given by Retarded Functional Differential Equations (RFDE) [13]. After applying certain simplifications, some special concepts arise such as models considering either discrete time delays, [14–16], distributed time delays, [1, 17, 18] or ordinary differential equations (ODE) of kinetic rates [2]. Although with varying performance, these models are shown to be capable of explaining the stability of biochemical networks [1, 13, 18–20].

Similar simplifications of RFDE as reported for parametric models can also be applied to hybrid semi-parametric models, i.e. either discrete delays or distributed delays of state variables in the kinetics or differential equations of the kinetics. The latter is not well suited for hybrid modelling, because neither kinetic function nor the kinetic values are known a priori and thus a solution or estimation of the kinetics is not straightforward. Distributed delays are also rather unlikely to be used, because one would have to introduce some function accounting for the delay, which is generally not known. Furthermore, some mathematical postulation of arbitrarily large delays for unknown weighting functions of the delayed variable would have to be assumed and this mathematical convenience is in limit biologically unrealistic (see [13]). Instead, the use of discrete delays in the inputs to the nonparametric structure is proposed herein. This is analogous to the application of discrete time series, namely **A**uto**R**egressive (e**X**ogenous), (AR(X)), models. This presents no limitation for application, as it is mathematically clear that a weighted discrete time series is equivalent to the integration of a time delay weighting function and thus analogous to the application of the distributed delay framework.

Unfortunately, the theoretically endless number of time lagged values of one variable would in practice lead to high computational times and identification problems of the network structure and parameters (see [21, 22]). In theory, an optimal number of time lagged values exists given by the ratio between redundancy and additional gain of information in the inputs. Several methods for identification of the optimal number of time lagged values such as autocorrelation, cross-correlation, or partial mutual information, have been proposed [21, 22]. However, they either assume that inputs are linearly correlated or are based on maximum information transfer and thus require known outputs, which unfortunately are not directly available in hybrid models.

Hence these methods cannot be applied here and thus the number of lagged values is rather chosen by trial and error, as done by several other authors [23, 24]. However, choice by trial and error is not a disadvantage, when (i) the delay is an important property of the system and when (ii) series of delays are systematically studied, since it can be expected that models that account for time delays perform the best when the studied and the "true" delays are congruent.

In this paper hybrid delay differential equations with discrete time series was determined to be a powerful method to identify delayed dynamics of ill-defined biochemical networks. This technique is described in detail in the results section. The technique was applied to a typical gene regulatory system where the transport of macromolecules between the cytosol and nucleus introduce strong delay dynamic effects. In addition, heterologous protein expression by recombinant *Pichia pastoris* was studied by assuming a hypothetical network with distributed time delays.