Field of the invention

[0001] The present invention concerns methods for calculating kinetic parameters of a reaction network which involves providing an intermediate objective function which comprises a linearized intermediate discrepancy function, wherein all hidden states have been eliminated, and minimizing the intermediate objective function, using a direct estimation method.

Description of related art

[0002] A chemical or biochemical reaction network describes the reactions, in general reversible, occurring between a set of reactants to form a set of products, wherein a reactant is a substance that takes part in and undergoes change during a reaction and a product is a substance that is the result of a reaction between one or more reactants. A reactant can simultaneously be a product and a product can simultaneously be a reactant in a reaction network.

[0003] The state of a reaction network is given by the set of concentrations of all reactants and products in the reaction network, often represented by a single concentration vector.

[0004] Each reaction in a reaction network can be described by a differential equation in the concentrations of the involved reactants and products, said differential equation comprising kinetic parameters that determine the reaction velocities among the reactants. The time evolution of the state of a reaction network is determined by the set of differential equations resulting from all reactions in the network. It should be understood that different reaction networks are described by different set of differential equations, comprising different kinetic parameters. [0005] Determining the kinetic parameters of a reaction network is of interest in chemical and biochemical applications. For example, in drug development the kinetic parameters of a reaction network comprising as reactants a drug and a drug target determine the suitability of the drug. In bioprocessing, the kinetic parameters of a reaction network are useful to optimize a bioreactor so as obtain a desired product with a minimal amount of the reactant substances and in a minimal amount of time.

[0006] The kinetic parameters of a reaction network are

determined from a sequence over time of measurements of the reaction network state, or a part of the reaction network state.

Wherein a measurement of a part of the reaction network state consists of a set of numerical values for one or more concentrations of the reactants and products in the reaction network state at a given time. Hereafter, a signal is a sequence over time of measurements of a part of the reaction network state, or of measurements of the reaction network state.

[0007] It is intuitively understood that a signal is acquired, or measured, by performing an acquisition method using a physical device (such as biochemical sensors). An acquisition method typically has physical limitations and the acquired signal contains noise and artefacts resulting from these physical limitations of the acquisition method.

[0008] A requirement on the signal is that the time interval between two consecutive measurements in a signal is smaller than the time scale of the reaction network for which the kinetic parameters are to be determined. The time scale of the reaction network is defined as the time the state of a reaction network needs to reach the equilibrium state, or the time of the shortest oscillation for the state of a periodic reaction network. The equilibrium state is the state in which the reaction network no longer changes over time and remains in the equilibrium state. The state of a periodic reaction network never reaches equilibrium but oscillates over time with a repeating pattern.

[0009] Acquisition methods performed to acquire a molecular interaction signal, molecular interaction which is described by a respective reaction network, rely on biochemical sensors. These biochemical sensors can be split into two categories; namely label-free sensors and sensors which operate using a marker. These sensors are generally known in the field of life sciences and are mainly used for the characterization of interactions of biologic or biochemical molecules. These characterizations commonly involve bringing two or more reactants, such as different types of sample molecules, into physical contact with each other for a set period of time and recording the signal for said period of time.

[0010] Acquisition methods using sensors which operate using a marker involve chemically attaching a marker, typically a fluorescent, absorbing or radioactive molecule, to a reactant or product which is to be detected. Here, the concentration of a state in a reaction network is determined indirectly by measuring the concentration of the marker attached to the reactant or product.

[0011] US 8,325,347 B2 describes an acquisition method which uses a label-free optical sensor to acquire a molecular interaction signal; the method is based on the waveguide interferometry where the local refractive index near the sensor surface is probed by an

evanescent wave of a waveguide mode and read out by

interferometry.

[0012] Acquisition methods which use label-free sensors include, but are not limited to, optical methods based on surface plasmon resonance (SPR) or waveguides, or methods based on surface- acoustic waves (SAW), thermal methods, or electro-chemical methods. Optical methods are based on the principle that biochemical molecules exhibit a different refractive index than an aqueous solution. Refractive index changes near the sensor surface result from the addition or

subtraction of molecules to the surfaces due to the interaction of molecules with either the sensor surface itself or another molecule attached to the surface. Using a resonant element- in case of SPR a metal layer supporting surface plasmons, or in case of waveguide sensors an optical waveguide supporting optical waveguide modes - the local refractive index changes can be then probed using an appropriate illumination and detection scheme, and the changes recorded in real time in order to measure the molecular reaction event. In this context, these changes correspond to the signal. Surface- Acoustic Wave and electro- chemical methods operate in a similar way, except that not the refractive index differences are physically

measured, but rather mass or permittivity differences.

[0013] An exemplary type of a known label-free optical sensor, which can be used in an acquisition method for molecular interactions, is a surface plasmon resonance-based biosensor (or SPR-based biosensor). SPR-based biosensors utilize a SPR based mass-sensing technique to provide "real-time" signals between a surface bound ligand and an analyte of interest. The SPR-based biosensor produces as output signal a sequence in time of density measurements in pg/mm2, density measurements which are proportional to the concentration of the product on the sensor surface.

[0014] In the abovementioned applications, and potential further applications, not all states in a reaction network can be measured. Hence, the reaction network state is split into two components namely the 'hidden state' and the 'observed state'. The 'hidden state' is a set of reactant and product concentrations in the reaction network that are cannot be measured due to physical limitations; and the 'observed state', is a set of reactant and product concentrations in the reaction network that can be measured.

[0015] A 'model signal' is the expected signal as determined by the differential equations of a reaction network model. [0016] Typically, the observed state is not observed directly but observed indirectly by a physically observable proxy state. The measured values over time for the physically observable proxy state are hereafter designated by observed signal. For example, the observed signal acquired by performing an acquisition method using a label-free biochemical sensor typically corresponds to a refractive index signal that is in (quasi-)linear relationship to the product concentration.

[0017] The observed signal is typically not equal to the model signal due to noise arising from the physical limitations of sensor(s) used in the acquisition method. The noise can be any combination of systematic deviations and stochastic deviations of the observed signal from the model signal.

[0018] A 'discrepancy function' is a mathematical expression for the discrepancies over time between an observed signal and a model signal or for the discrepancies over time of an observed signal, or its integral or derivative, with respect to a set of differential equations.

[0019] An 'objective function' is a function that summarizes a discrepancy function over time. The said objective function comprises at least a discrepancy function as parameter, where said discrepancy function comprises kinetic parameters. The objective function can further comprise one or more additional parameters. The values of the kinetic parameters that minimize said objective function are the values that best explain an observed signal in terms of a reaction network.

[0020] The procedure of finding the minima of the objective function is called minimization, where said minimization is performed using a mathematical optimization method.

[0021] Any minimization problem of an objective function is equivalent to the maximization of the opposite of the same objective function, mathematically expressed as

In this application the term minimization will be used to refer to minimizing the discrepancies of an observed signal with respect to a reaction network model; however, it should be understood that every minimization step mentioned in this application could be replaced with a mathematically equivalent maximization step (since the minimization can always be rewritten as a mathematically equivalent maximization).

[0022] State of the art computation of kinetic parameters from an observed signal uses iterative fitting methods to minimize an objective function (Jesudason, 201 1 ; Myszka, He, Dembo, Morton, & Goldstein, 1998; Voit, Martens, & Omholt, 2015). For example, (Myszka & Morton, 1998) use the Levenberg-Marquardt nonlinear minimization

algorithm, and (Canedo & Gonzalez-Hernandez, 201 1) use a second order gradient optimization algorithm. When the underlying reaction network model is not analytically solvable (e.g. mass transport model) an additional numerical integration scheme is needed to solve the differential equation(s). The repeated solving of the reaction network model, including necessary numerical integrations, is computationally intensive and delays the workflow of evaluating an observed signal. In particular, the workflow becomes inconveniently long when repeated evaluation is needed to tune the fit. Examples of parameters to tune are starting times and end times determining the observed signal range used for the computation. Further on, computing the

estimation errors of the obtained values for the kinetic parameters is a complex problem (Johnson, Simpson, & Blom, 2009), which needs repeated fitting and becomes unfeasible on a single processor.

[0023] The simplest reaction network is the Langmuir model with the reversible reaction

A + B ¹ AB, (2·) where an analyte A and an immobilized ligand B reach equilibrium with the product AB. The reactants A and B associate to form the product AB, while simultaneously the reactant AB dissociates to form the products A and B. Hence, A, B and AB are simultaneously reactants and products.

[0024] Mathematically, the reaction network is described by the differential equations

where [A], [B] and [AB] are the concentrations of A, B and AB. The velocity of [A] and [B] associating to form the product [AB] is determined by the kinetic parameter k _a , and the velocity of [AB] dissociating to form the products [A] and [B] is determined by the kinetic parameter k _d .

[0025] An acquisition method is then performed to obtain an observed signal which, indirectly, represents the product

concentration [AB]. When performing an acquisition method using a label-free sensor, the product concentration [AB] is observed indirectly via a refraction index signal. Importantly the analyte concentration [A] is held constant when performing the mentioned acquisition (the analyte concentration [A] is held constant for the duration of the refraction index signal). Hence, [A] and [AB] are herein observed states. In contrast, the ligand concentration [B] cannot be observed and is herein a hidden state. A sequence of different analyte

concentrations results in a sequence of signals. In this example a label- free acquisition method is performed.

[0026] Despite the hidden state, the kinetic parameters are entirely determined by the two observed states. To eliminate the hidden state [B], the differential equations (2. a) and (2.b) are added and their sum is then integrated to obtain d[AB] d[B]

= 0 = [B] = R _max - [AB],

dt dt max where the integration introduces the 'integration constant' R _max (or equivalently 'constant of integration'). This 'integration constant' is an additional non-kinetic parameter that needs to be determined.

[0027] Substituting the result of Equation (4.) in Equation (3.a) yields the differential equation which is a differential equation in the observed product concentration [AB] only.

[0028] During the association phase, the analyte has a constant concentration [A] = C, and the solution to the differential equation (5.) is given by where R(t) is the model signal of the Langmuir model during the association phase.

[0029] During the dissociation phase, the analyte concentration is [A] = 0, and the solution of the differential equation is given by

R(t) = AB = R ₀ exp (— fc _d t) , (7·) where R(t) is the model signal of the Langmuir model during the dissociation phase, and R ₀ is the model signal at the end of the association phase. [0030] Figure 1 show baseline observed signals, association observed signals and dissociation observed signals for the Langmuir model at three different analyte concentrations and a fixed initial ligand concentration. Typical extensions to the Langmuir model are the mass transport, the heterogeneous ligand, the bivalent, the conformational change, and the heterogeneous analyte models.

[0031] The aim is to calculate kinetic parameters of a reaction network model so as to best explain the observed signal R _0bserved (f) with the model signal R _m0dei J _> V) which in state-of-the-art

applications is achieved by using the discrepancy function

The model signal R _m odei depends on the parameters p that comprise the kinetic parameters.

[0032] State-of-the-art applications typically use for objective function the chi-squared expression which summarizes the discrepancies over time into a single scalar value obtained by summing over time the square of the discrepancy function values. [0033] The model signal P _modei (t, p) of the Langmuir model during association is given by Equation (6.), and during dissociation is given by Equation (7.). The parameters of the Langmuir reaction network are p = (k _a , k _d , R _max ), where k _a and k _d are the kinetic parameters.

[0034] Methods used in state-of-the-art applications to minimize the objective function are of iterative type, including gradient based methods such as the Gauss-Newton algorithm, the Levenberg- Marquardt algorithm, and variants or extensions thereof.

[0035] The major drawback of iterative minimization methods, in particular when the model requires numerical integration, is the significant computation time required to achieve convergence. The computation time introduces significant delays in the experimenters' workflow who is unable to quickly evaluate the kinetic parameters for different configurations.

[0036] Further on, the computational complexity of iterative methods prevents real time Monte-Carlo simulations to compute the estimation error, tune hyperparameters such as the fit range, or perform model selection.

[0037] An alternative to iterative methods for determining the kinetic parameters are direct estimation methods that minimize an objective function summarizing the discrepancy of an observed signal with respect to a different discrepancy function. The discrepancy function used in direct estimation methods is defined by the

differential equations of the reaction network. Direct estimation methods are methods that obtain the desired result in a single closed form computation and do not require iterative steps that converge to the solution without reaching it exactly. A common direct estimation method is ordinary least square estimation to determine the variables (i.e. desired result) for a linear system of equations.

[0038] The direct estimation method derived by Kovacs & Toth (2007) is given by for the reaction network equation x(t) = F(x(t)) · k + F _r (x(t)), x(0) = x ₀ , (11.) wherein the variable x(t) is a vector, and the differential equation implicitly defines multiple differential equations written in vector form.

[0039] The Equations (3.a) and (3.b) are a special case of the reaction network Equation (11.), which becomes apparent when writing Equations (3.) in the form

[0040] The kinetic parameters k _a and k _d can then be determined directly by computing the direct estimator in Equation (10.), estimator which does not require any iteration to converge to the desired result.

[0041] When the model signal equals the observed signal, and the observed signal is sampled at a time resolution smaller than the rate of change, the difference between the kinetic parameters determined with an iterative method and a direct estimation are negligible for practical purposes. When discrepancies between the observed signal and the model signal arise, significant differences can occur between the kinetic parameters determined with an iterative method and a direct estimation method. The differences depend on the statistical properties of the discrepancy values, and these differences are not mathematically understood in the general case.

[0042] The requirement of Equation (10.) to observe simultaneously the product concentration [AB] and the ligand concentration [B] rules out this direct estimation method as an alternative to determine kinetic parameters from observed signals where only the product concentration is observed, and the ligand concentration is a hidden state of the reaction network.

[0043] A solution to the limitation of direct estimation methods to signals observing all states in a reaction network has recently been proposed (Dattner, 2015). The solution involves determining an adequate functional basis to describe the observed states and hidden states. For example, in Equations (12.) and (13.) the hidden ligand state can be approximated by

where the (<p ₀ (t), ... , f _k (ΐ)} form a finite functional basis that allows for a good approximation of the ligand concentration over time (e.g. polynomial basis { ί, c, c ² , ... }), and /? _fc are parameters to be determined. In vector form, any state can be written as state (t) = fi _state f(ΐ), (15.) the parameter vector and <p(t) = tional basis vector.

[0044] For observed states, the parameters _state can be

determined by a regularized linear regression estimator state(t), (16.) wherein l is a regularization constant; the value of l is pre-defined by a user; t is the vector of times at which the observed states are observed and f( ) is the matrix with columns given by finite

functional basis vectors at the times t, f(ΐ) is the transposed matrix of f( ) and the matrix inverse of the matrix product f(ΐ)f(ΐ) ^T . [0045] For the example, the kinetic parameters can then be determined as a function of the parameters b _B and b _AB , which relate to the states using the direct estimation method of Equation (10.) as

The parameters b _AB are determined from the observed state [AB](t) using Equation (16.). The parameters b _B cannot be determined because the corresponding state [S](t) is a hidden state. However, initial values for b _B must be provided for the subsequent minimization of the objective function in Equation (18.). Initial values for [S](t) are obtained from the relation [B](t) = [B]( 0) - [AB](t), wherein [B]( 0) is an initial ligand concentration and wherein a value for [£]( 0) is estimated based on the maximum value over time of the observed state [AB](t). The initial values of the parameters b _B are determined by applying Equation (16.) to the initial values for the state [S](t) = [b](0) - [AB](t) wherein [£]( 0) is replaced by the estimated value.

[0046] Solving for the kinetic parameters then becomes a minimization problem with respect to the parameter matrix /? _state , where the matrix element /¾ _ate corresponds to the weight of the yth basis function for the ith state. The parameter matrix can be written in the bi-partite form b _3ίaί6 = \bh _ίaabh> b o _b served where the hidden and observed states are made explicit.

[0047] For the example, the parameter matrix is given by /? _state =

[bB. bABΐ

[0048] The objective function of the minimization problem reads

where the first term enforces a solution to the reaction network differential Equation (11.) and the second term minimizes the chi- squared between the solution and the observed states. The parameter ^ _di // controls the relative weighting of the differential equation discrepancy and the discrepancy between the solution and the observed state. This parameter is initialized with = 0.1 and increased each time a solution to /? _state is found, until the value of be _ί a _ίb ^no longer changes when increasing

[0049] This minimization problem has the advantage of being dominantly linear in the parameters /? _state and not to require integration of the differential equations, which is computationally much more efficient to solve. However, a major drawback of the method is that it remains iterative and the determined kinetic parameters are only an approximation whose precision depends on how well the functional basis approximates the true solution of the differential equation.

[0050] Besides the limitation of direct estimation methods to observed signals that observe all states in a reaction network, the statistical properties of current direct estimation methods are not always optimal. An undesirable statistical feature that can arise in iterative methods and direct estimation methods are biases in the determined kinetic parameters. A bias arises when the expected value of a result differs from the true underlying quantitative value being estimated. Biases are difficult to correct because they require knowledge of the distribution underlying discrepancies between the observed signal and the model signal. For Gaussian independent noise in the observed signal, direct estimation methods are known to be more prone to biases and have lower statistical efficiency than iterative methods. The statistical efficiency being the number of observations needed to achieve a given performance in computing the kinetic parameters.

[0051] However, recent mathematical progress has led to the development of a bias corrected direct estimation method (Holte, 2016), robust to Gaussian time independent noise in the observed signal. The asymptotic consistency of such corrected direct estimation methods has been proven for independent observation errors (Dattner & Gugushvili, 2015; Vujacic, Dattner, Gonzalez, & Wit, 2015).

Nonetheless, the statistical properties of direct estimation methods remain unknown for dependent observation errors and heavy-tailed observation errors, which are common in applications.

[0052] In summary, state-of-the-art iterative methods which are used to determine kinetic parameters of a reaction network are insufficient because they are computationally intensive and delay the workflow of evaluating observed signals. State-of-the-art direct estimation methods which are used to determine kinetic parameters of a reaction network are insufficient because they cannot be applied to observed signals that do not observe all states of a reaction network. The recent development of a direct estimation methods that applies to reaction networks with hidden states (Dattner, 2015) relies on a parametric functional approximation of the hidden states that is optimized using iterative methods; however disadvantageously this approach results only in an approximation and requires case by case educated guesses for the choice of the functional basis.

[0053] It is an aim of the present invention to obviate, or mitigate, at least some of the above-mentioned disadvantages which are associated with existing methods for calculating kinetic parameters of a reaction network.

Brief summary of the invention

[0054] According to the invention, these aims are achieved by means of a method having the features recited in the independent claims; wherein the dependent claims recite optional features of preferred embodiments. [0055] Advantageously, the method of the present invention can be applied to observed signals that do not observe all states of a reaction network; additionally, the method of the present invention is faster than conventional iterative methods, enabling advanced statistical analysis in real time, such as estimating errors bounds in the parameters, selecting the data range with relevant information about the kinetic parameters of the reaction network, and performing model selection.

[0056] Advantageously, the method of the present invention can be used to efficiently calculate kinetic parameters of a reaction network from an observed signal or sequence of observed signals of said reaction network.

[0057] It should be understood that any definitions of terms provided in the 'Description of related art' section also apply to the present description.

Brief Description of the Drawings

[0058] Figure 1 shows baseline observed signals, association observed signals and dissociation observed signals, acquired using an SPR-based biosensor, for a product concentration in the Langmuir reaction network at different analyte concentrations and a fixed initial ligand concentration.

[0059] Figure 2 shows the observed signals from Figure 1 with overlaid association model signals and dissociation model signals obtained with the kinetic parameters calculated with an embodiment of the present invention.

Detailed Description of possible embodiments of the Invention [0060] The invention will be better understood with the aid of the description of the following embodiments which are given by way of example only.

[0061] A reaction network (such as a reaction network arising in chemistry, biochemistry, pharmacokinetics, and related fields) is described by a set of differential equations; in this example the reaction network is described by a set of first order differential equations each differential equation of the set being of the form where x(t) = (x·^ t), ... , x _w (t)} e ^N is the reaction network state, t is time, Si _j e E are the kinetic parameters of interest, and r _ijk are stochiometric coefficients of the reaction network. The values of the stochiometric coefficients are typically determined by the law of mass action. It should be understood that different reaction networks are described by different sets of differential equations.

[0062] Typical acquisition methods output a signal that observes only some of the concentrations in the reaction network state x(t) over time t. A state (e.g. the concentration of a component in the reaction network) in the reaction network is observed by means of a physically observable proxy state (e.g. a refractive index) which are representative of that state, such a state is an observed state; for other states, no physically observable proxy state which are representative of that state are available, such a state is a hidden state. Hereafter, to distinguish between the observed states and the hidden states in the reaction network, the reaction network state is written in the bi-part form x(t) = (x ₀ (t), X _f t(t)}, which is the concatenation of the observed states x ₀ (t) and the hidden states x _h (f).

[0063] In the present embodiment, in order to avoid the costly computation of hidden states, the hidden states are eliminated from said set of differential equations, by substituting parameter(s) which represent hidden states with equivalent expression(s) comprising parameter(s) representing observed states, so as to form one or more intermediate differential equation(s) which are without hidden states.

[0064] To obtain said one or more intermediate differential equation(s), said set of differential equations is augmented with derivatives or integrals of said differential equations. The elimination of the hidden states may introduce higher order derivatives of the observed state into the one or more intermediate differential equations. The higher order derivatives can include second order derivatives of the observed states and/or third order derivatives of the observed states. The one or more intermediate differential equations depend only on the observed states x ₀ (t).

[0065] The one or more intermediate differential equations can admit a vector space of solutions larger than the space of solutions to the original kinetic system of differential equations. For example, the intermediate differential equation can have solutions with negative response, which is not physically possible. To obtain the unique solution that describes the reaction network, a set of ri _j initial conditions {f _itial (p, e [), - ,d _tial (p, e _n .)}, specific to the reaction network, are preferably enforced. Where the parameters {e[, , e _n .} are slack variables which accommodate discrepancies of the

observations from the model initial conditions. The slack variables are to be minimized.

[0066] In this example the one or more intermediate differential equations define a discrepancy function f _discrepancy (J v)·

[0067] The discrepancy function fdiscrepancy (J _> V ) comprises parameters p = {p ₁ ... , p _n written in vector form, said parameters p comprising the kinetic parameters of said reaction network; the parameters p of the discrepancy function may further comprise

'constants of integration' resulting from the elimination of hidden states from said set of differential equations. The expression 'constant of integration' is a terminology commonly used in the art, but herein, the 'constants of integration' are unknown parameters (i.e. they are not provided constants), which are not kinetic parameters; in this example a value for any 'constants of integration' present in the discrepancy function is determined in a subsequent step.

[0068] An objective function is then defined. The discrepancy function fdiscrepancy (J _> V ) is ^a parameter of the defined objective function. In this example, a penalized chi-squared expression

(Equation 20.) defines said objective function: wherein fdiscrepancy (J _> V ) is the discrepancy function, p are the parameters of the discrepancy function, and L = diagC^, ... , A _n ) is a diagonal penalty matrix. p ^T A p defines an additive penalty term which penalizes each parameter p _t with a respective penalty coefficient L _έ . The penalty coefficients are additional parameters. In this example, the discrepancy function / _discrepancy (t, p) and additive penalty term p ^T A p, define the parameters of the objective function. However, it should be understood that the additive penalty term p ^T A p is not essential to the invention; the objective function may have the discrepancy function / _discrepancy (t, p) only as a parameter.

[0069] The discrepancy function f _discrepancy ( _. t, p ) is non-linear in the parameters p, which prevents the use of a direct estimation method to determine the values of the parameters that minimize the objective function. In order to determine the values for the parameters that minimize the objective function, the objective function is

reparametrized so as to linearize the discrepancy function.

[0070] The discrepancy function f discrepancy ( _. t, p ) is linearized by reparametrizing each of its parameters p with a bijective reparameterization function /^ _parameteri z _ation that maps expressions in the parameters p = {p ₁ ... , p _ni } to intermediate parameters k =

{ki, ..., k _n2 } as

The number of intermediate parameters is equal to, or greater than, the number of parameters in the discrepancy function (i.e. n ₂ > n^).

The inverse of the reparameterization function defines expressions in the intermediate parameters k that map to the parameters p as

P = / reparameterization (k). (22.)

The reparameterization of the parameters p of the discrepancy function f _d iscrepancy (J _> V ) ' ^s done by mapping, using the inverse of the reparameterization function / _p a _rameteri z _ation ^eac h of the parameters in p in the f _d iscrepancy (J _> V ) t° ^a respective expression of intermediate parameter in k ; and then substituting each of the parameters p in the its respective expression in intermediate parameters in k. It should be understood that discrepancy functions of different reaction networks are linearized by different

reparameterization functions.

[0071] Once each of the parameters p of the discrepancy function have been substituted with their respective expression of intermediate parameters in k, this yields a linearized intermediate discrepancy function that is linear in the intermediate parameter k. In other words, the reparametrized discrepancy function defines a linearized intermediate discrepancy function.

[0072] In this example, the objective function further had the additive penalty term p ^T A p as a parameter. The additive penalty term p ^T A p is also reparametrized, by mapping each of the parameters in the additive penalty term to its respective expression in

intermediate parameters using the inverse reparameterization function / _{reparameterization} ϊ and then substituting each of the parameters in the additive penalty term with its respective expression in

intermediate parameters. As mentioned, the additive penalty term p ^T A p is not essential to the invention, accordingly the step of reparametrizing the additive penalty term p ^T A p is not essential to the invention.

[0073] Thus, as described above, the present invention involves a step of reparametrizing the objective function; specifically, if the objective function has only the discrepancy function f _di scre _p anc _y (.t, p ) as a parameter then the step of reparametrizing the objective function involves reparametrizing the discrepancy function f _di scre _p anc _y (.t, p ) only; if on the other hand the objective function has the discrepancy function f _di scre _p anc _y (J _> V ) ^{as a} parameter and has one or more other parameters (such as the additive penalty term p ^T A p mentioned in the above example) then the step of reparametrizing the objective function involves reparametrizing the discrepancy function

and reparametrizing said one or more other

parameters.

[0074] Accordingly, in the above example, the objective function, now having a linearized intermediate discrepancy function (i.e. the reparametrized discrepancy function f _di scre _p anc _y (.t, p )) as a parameter, and the reparametrized additive penalty term, defines an

intermediate objective function.

[0075] The linearized intermediate discrepancy function is a vector with components wherein the index i indexes the i components of the vector, and where each component is linear in the intermediate parameters k = are each observation functions of observed states in the reaction network over time and derivatives of the observed state over time, wherein m is the index of the highest order derivative in said intermediate differential equations. The observation functions result from the elimination of the hidden state from said set of differential equations. It should be noted that the range of the indices i and j in the linearized

intermediate discrepancy function (i.e. intermediate differential equations) may be different from the range of the indices in the set of differential equations.

The linearized intermediate differential equations that define the linearized intermediate discrepancy function consists of one or more differential equations. Each of the one or more differential equations corresponds to a respective component of the vector representing the linearized intermediate discrepancy function, said components being defined by Equation (23.).

[0076] {X _O (T _L ), ... , x ₀ (t _T )} is an observed signal, or sequence of observed signals, which represents an observed state within the reaction network at times {t ₁ ..., t _r ); the values over time of the observation functions (t)) are entirely determined by values of the observed signal and the higher order derivatives of the observed signal.

[0077] In a preferred embodiment, an acquisition method is performed so as to obtain an observed signal which directly or indirectly represents an observed state, or a sequence of observed signals which directly or indirectly represent a sequence of an observed state. Preferably a label-free acquisition method is performed; and more preferably the acquisition method is performed using optical sensors with high sensitivity, large measuring range and high readout speed, particularly consisting of integrated-optical waveguides and a readout device, as they find use, for example, in pharmacology or in diagnostics.

[0078] For example, a label-free acquisition method may be performed, by attaching (immobilizing) of one or more "ligands" (a reactant such as antibodies or drug targets) to a solid support on a sensor surface (i.e. a sensitized surface on a sensor chip which is adapted to be read out by a detection scheme which outputs the measurements over time that compose the recorded observed signal). The fluid conduit containing the sensor surface is generally called flow cell or fluidic chamber and allows bringing a fluid containing the other reactant(s) to be investigated (analyte) to be brought into contact with the ligand. Thereby molecule(s) to be investigated

(analyte) have the opportunity to react with the immobilized ligands at the solid support on a sensor surface and an eventual product concentration is measured and characterized. A typical molecular interaction signal is acquired by first contacting the surface-bound ligand first with a neutral buffer solution in order to establish a base signal without reaction ("baseline"), followed by contacting the surface-bound ligand with a fluid containing the actual analyte or sample (such as an antigen) so that the surface-bound ligand reacts with the analyte or sample, and the association phase where a reaction occurs can be monitored, and optionally followed by contacting the surface-bound ligand again with a neutral buffer solution in order to monitor the dissociation phase of the analyte or sample by removing the analyte from the flow cell. In other words, the concentration of the analyte is increased in a step-wise fashion for characterizing the association phase and decreased again in a step wise fashion for characterizing the dissociation phase. Typically, the aforementioned steps are repeated for several concentrations of analyte. Typically, during the whole time of the experiment, measurements are recorded, resulting in an observed signal or sequence of observed signals, which may be further analyzed with the present invention to obtain the kinetic parameters.

[0079] An exemplary sequence of observed signals, which represent an observed state that has been obtained by performing an

acquisition mention using a SPR-based biosensor, is shown in Figure 1. In this example the observed state represented by the signal is the concentration of the product resulting from an analyte reacting with a ligand. The y-axis indicates the measurement value (here in pg/mm2) and the x-axis indicates the time (here in seconds). Initially, buffer is passed over the measurement surface giving the "baseline signal". During sample injection, an increase in the observed signal is due to the reaction between the ligand and the analyte until the observed signal plateaus when the equilibrium state is reached. The

measurements acquired during sample injection give the "association signal". At the end of sample injection, the sample is replaced with a continuous flow of buffer and a decrease in the observed signal reflects the dissociation, or release, of analyte from the surface. The measurements acquired during the continuous flow of buffer give the "dissociation signal". The slope of the association/dissociation curves provides valuable information regarding the kinetic parameters, and the observed signal represents surface concentration (i.e., the concentration of the product is related to the change in product density on the surface). The observed signal(s) obtained by performing an acquisition method using biosensor(s) based on other detection principles will have a similar appearance.

[0080] Typically, a reaction network is characterized by a sequence of observed signals. A typical sequence of observed signals is characterized by repeatedly acquiring a baseline signal, an association signal, and a dissociation signal for different concentrations of a same analyte. Wherein any or several of the steps of acquiring a baseline signal, an association signal, or a dissociation signal can be skipped. The lowest and highest measured analyte concentrations define an analyte concentration range. The kinetic parameters only depend on the analyte and ligand substance but do not depend on the analyte concentration. The aim is to determine the kinetic parameters values of a reaction network that best explain the entire sequence of observed signals, said observed signals being acquired at different analyte concentrations. Determining the said values of the kinetic parameters is therefore done using the entire sequence of observed signals, or a suitable sub-sequence of the observed signals. Hereafter, an observed signal is used interchangeably to designate a single observed signal or a sequence of observed signals.

[0081] In another embodiment, the observed signal or sequence of observed signals is obtained performing other acquisition methods to acquire a molecular interaction signal, such as acquisition methods with sensor(s) that operate using a marker.

[0082] In another embodiment, the observed signal or sequence of observed signals is obtained with an acquisition method used in the field of bioprocessing, which are important acquisition methods in a wide variety of industries such as pharmaceutical, genetics, food, ecology and water treatment. In bioprocessing acquisition methods (and techniques) a measurement of the reaction network state, or part of the reaction network state, relies on optical chemical sensor technologies, as for example the optical chemical sensor technology described in (US 7,041,493 B2). In such optical chemical sensor technologies, an excitation source produces light which excites an optical chemical sensor to generate emission and/or cause absorption. The emission and/or absorption is measured by a detector. The luminescence emitted from the chemical sensor or the amount of light absorbed by the chemical sensor is related to the concentration of an analyte. Such as oxygen. If the luminescence emitted changes, or if the amount of light absorbed changes, then the concentration of the analyte has changed. Further examples of measurements are carbon dioxide concentration, biomass concentration, oxygen concentration, substrate concentration, and glucose concentration.

[0083] The one or more observed signal(s) (representing observed states within the reaction network), obtained with an acquisition method, are used to form an observation matrix as follows:

[0084] For given observed signal(s) over time (x ₀ (ti), ... , x ₀ (t _x )}, the values of the observation functions (t)) over time are summarized in an observation matrix, as shown in Equation (24.) below.

Each column in the observation matrix corresponds to a respective intermediate parameter. The first column corresponds to the intermediate parameter k _lt the second column to the intermediate parameter k ₂ , and so on. The rows in the observation matrix

summarize the components of the intermediate discrepancy function and the individual measurements of the observed signal. The first consecutive i rows (where i is the number of components in the intermediate discrepancy function) correspond to the i components of the intermediate discrepancy function at time t _lf the following i consecutive rows correspond to the i components of the intermediate discrepancy function at time t ₂ , and so on.

It is to be understood that an equivalent observation matrix and observation vector can be obtained by: permuting the rows are equivalent; or permuting the columns of the observation matrix and the respective intermediate parameters.

Analogous to the observation matrix, the values of the observation function £ _j (x ₀ (t), ... , d _m x ₀ (t)) are summarized in an observation vector b.

[0085] The observation matrix and observation vector are typically built from a sequence of observed signals. For example, reaction networks representing a molecular interaction signal have for observed state an analyte [A] and a product concentration [ AB ]. The analyte concentration [A] is given a sequence of provided values (i.e. provided by the user) each imposed to the reaction network for a time interval T _t . A positive analyte concentration interval give rise to an association observed signal and an interval with an analyte concentration of zero gives rise to a dissociation observed signal. The sequence of analyte

concentrations gives rises to a respective sequence of observed signals.

[0086] To build the observation matrix and observation vector from a sequence of observed signals, the first step comprises determining for each observed signal, which are indexed by i in the sequence of observed signals, the respective interval observation matrix X _t and interval observation vector b _t .

[0087] The observation matrix of a single observed signal is composed of two parts: a shared interval observation matrix X _{i 0} related to the intermediate parameters k ₀ shared by all observed signals in the sequence, and an interval specific observation matrix X _{t l} related to the interval intermediate parameters k _{t l} that are specific to that interval.

[0088] Interval specific intermediate parameters occur for example in observed signals obtained by performing a label-free acquisition method that use a refractive index as observable proxy signal. A refractive index signal is often susceptible to bulk refractive index mismatches, such as small differences in refractive index between the neutral buffer solution and the fluid containing the sample. Such mismatches can cause an observed signal offset, which is an offset constant in time by which the observed signal differs from the true underlying observed state during that time interval. Each interval defines a single observed signal. To accurately compute the kinetic parameters from such an observed signal, the observed signal offset needs to be computed jointly with the kinetic parameters. Therefore, an interval intermediate parameter that represent an observed signal offset is introduced for each interval

[0089] The observation matrix, the observation vector, and the intermediate parameters for a sequence of observed signals are given by

where X _{i 0} is the shared interval observation matrix of interval i, X _{t l} is the interval specific observation matrix of interval i, b _t is the

observation vector of interval i, k ₀ are the shared intermediate parameters and k _{t l} are the interval specific parameter(s) of interval i.

It should be noted that the case of determining kinetic parameters from a single observed signal is done by using a sequence of observed signals of length one, where the single observed signal is the only observed signal in the sequence.

[0090] The observation matrix and observation vector are used to form a compact representation of the components of the intermediate discrepancy function over time, as given by Equation (27.) intermediate discrepancy function over time = Xk + b (27.)

Xk + b defined a compact representation of the intermediate discrepancy function; advantageously this compact representation of the intermediate discrepancy function can be used to simplify subsequent steps.

[0091] Next, numerical values of said observation matrix and said observation vector are determined. In this example, in order to determine the numerical values of said observation matrix and said observation vector, the derivatives (dx ₀ (t), ... , d _m x ₀ (t)} of the observed signal which represents an observed state within the reaction network are computed. The derivatives of the observed signal can be estimated with methods known in the art, such as finite differentiation, local polynomial regression, convolutional filters, or a kernel smoother.

[0092] In one embodiment of the present invention, the derivatives of the observed signal are computed using the central finite difference scheme

[0093] In a preferred embodiment, a smoothed observed signal and smoothed observed signal derivative(s) are used instead of the observed signal, respectively instead of the observed signal

derivative(s). The smoothing alleviates the drawback of finite difference schemes, including the central finite difference scheme, which is their sensitivity to noise. The smoothing of the observed signal and observed signal derivative(s) can be performed by methods known to a person skilled in the art, such as local polynomial regression methods, kernel smoothing methods, smoothing spline methods, or convolutional filter methods.

[0094] In an embodiment of the present invention, the smoothed observed signal and derivatives of the smoothed observed signal are estimated using smoothing splines. The smoothing spline method estimates, in a first step, a smoothed observed signal satisfying the equation where f(t ) is the value of the smoothed observed signal at time t and l is a penalty term controlling the smoothness. The values of the smoothed signal can be computed from the observed signals in closed form using methods known in the art. In a second step, the

interpolating splines are then computed from the smoothed observed signal. The derivatives of the smoothed observed signal are then estimated at a given time by computing the derivatives of the interpolating spline at that time. In a subsequent direct estimation, the smoothed observed signal is used in place of the observed signal and the smoothed derivative estimates are used in place of the derivative estimates obtained in Equation (28.).

[0095] Next the compact representation of the linearized

intermediate discrepancy function Xk + b is substituted for the linearized intermediate discrepancy function which is a parameter of the intermediate objective function, so that the intermediate objective function becomes a compact intermediate objective function:

® 2 ® ®

Xintermediate ^~ (Xk + ft) + k A k, (30.) wherein X is the observation matrix, ft is the observation vector, k are the intermediate parameters in vector form of the linearized

intermediate discrepancy function, and A' is a diagonal matrix of intermediate additional parameters that penalize the intermediate parameters k of the linearized intermediate discrepancy function.

Equation (30.) defines a compact intermediate objective function.

[0096] The minimum with respect to the intermediate parameters of the compact intermediate objective function defined in Equation (30.) is found by setting the gradient of the compact intermediate objective function with respect to the intermediate parameters to zero, as given by the equation wherein X is the observation matrix, b is observation vector, A' is the diagonal matrix of intermediate additional penalty parameters, and k are the intermediate parameters.

[0097] The values of the intermediate parameters k that satisfy Equation (31.) are given by the direct estimation method which are the values that minimize the compact intermediate objective function. The values of the diagonal matrix of intermediate additional penalty parameters (L') in the direct estimation method of Equation (32.) are pre-defined by the user. In a simple embodiment, the user can set the diagonal matrix of intermediate additional penalty parameters to zero. In further embodiments, the user can determine the diagonal matrix of intermediate additional penalty parameters using methods known in the art, such as the Tikhonov regularization method, or the cross-validation method.

[0098] In the present example, a direct estimation method is used to determine values of the intermediate parameters k that satisfy Equation (31.); it should be understood that any suitable direct estimation method could be used. Preferably a robust least square regression estimation method is used.

[0099] The observed signal typically has a variety of quality issues as for example baseline slopes, spikes from air bubbles, oscillations and jumps originating from the physically observable proxy signal obtained by performing an acquisition method using biosensor(s) or sensor(s) used in bioprocessing. Such systematic deviations of the observed (proxy) signal from the model signal introduce strong auto correlations and heteroscedasticity in the residuals. To mitigate the impact of quality issues in the observed signal, further embodiments determine values of the intermediate parameters k that satisfy

Equation (31.) with a robust and bias corrected direct estimation method; this minimizes statistical biases and estimation errors in the estimated intermediate parameters, said biases and estimation errors depending on the statistical structure of the observed signal. In these further embodiments the direct estimation method which is used to determine values of the intermediate parameters k that satisfy

Equation (31.) can be any other known direct estimation method, such as the generalized least square estimation method, the iteratively reweighted least square estimation method, least absolute deviation estimation method, Bayesian linear regression estimation method, and other variants of robust closed form estimation methods.

[00100] In a further embodiment, the parameters are constrained to a physically plausible range with a provided lower and upper bound. Hard constraints on the parameters are implemented as follows. First, the constraints on the intermediate parameters are determined by applying the reparameterization function to the bounds of the range to which the parameters are constrained. Then the constraints on the intermediate parameters are enforced by using a direct estimation method such as non-negative least square regression or bounded- variable least squares, methods known to a person skilled in the art. The enforcement of the constraints on the intermediate parameters implies that the constraints on the parameters are satisfied as well.

[00101] As mentioned, the values of the intermediate parameters k that satisfy Equation (31.) are the values that minimize the compact intermediate objective function; and these values are determined using a suitable direct estimation method. Once the values for the intermediate parameters /r that minimize the compact intermediate objective function have been determined, the values of the

parameters p that minimize the objective function are then

determined.

[00102] The values of the parameters p that minimize the objective function are determined by computing for each parameter p the value of the respective expression in the intermediate parameters k, said expression being determined from the inverse of the

reparameterization function /p _{arameterization} -

[00103] When the intermediate differential equation admits a vector space of solutions, an additional set of initial conditions

{ filitiai ip ) = 0, - ,/¾ _tial ( ) = 0} are preferably enforced, because the values of the parameters p cannot be uniquely determined from the inverse of the reparameterization function p _{arameterization} - Instead, the values of the parameters p are determined by the minimization problem argmin[/ _r reparameterization (p) - k where p are the parameters to be determined, k are the intermediate parameters with known values, [e[, , e .} are slack variables to accommodate discrepancies of the observations from the model initial conditions, {/P _iai (pp), - ,fi _iaI (p, _£ )} is the set of n _t initial conditions for the parameters, and {P ··· ,X _n. } are penalty parameters for the slack variables, and values for {P ··· ,X _n. } are pre-defined by a user. The minimization can be solved using standard optimization methods such as Levenberg-Marquardt. In further embodiments, the user can determine values for the penalty parameters using methods known in the art, such as the Tikhonov regularization method, or the cross-validation method. This optimization has a very low

computational cost compared to fitting the complete model to all data points.

[00104] In said step of determining the parameters p from the intermediate parameters k, applying the inverse of the

reparameterization function /p _{arameterization} t° the intermediate parameters k may introduce additional biases to the parameters p. The present invention may optionally, include a step of correcting biases in the parameters p ; correcting biases in the parameters p may be done using standard expectation computation that relies on an assumption for the probability distribution of the intermediate parameters k. The probability distribution of the intermediate parameters may be derived from the distribution of the values of the intermediate discrepancy function over time using known methods such as an analysis of variance.

[00105] Assuming that the intermediate parameters k follow a probability distribution p(k), the expected value of the parameters p is given by

E[P] = J¾ fr reparameterization ( k ) p{k)dk. (34.)

[00106] The bias corrected value of the parameters p is given by

Pcorrected /reparameterization (k estimated ) + /reparameterization (E[k]) - E[p], (35.) wherein k _estimated is the value, and E\k] the bias corrected value, of the intermediate parameters k and E[p] is the expected value of the parameters for the chosen distributional assumption.

[00107] In an embodiment where the observed state has

discrepancies with respect to the estimated reaction network differential equation, the estimated kinetic parameters do not coincide with the kinetic parameters corresponding to the solution of the reaction network model differential equations that minimizes the chi-squared with respect to the observed state. To obtain the kinetic parameter corresponding to the least square solution of the reaction network model with respect to the observed state, the estimated parameters are preferably further corrected using known methods in the art (e.g. the method by (Dattner, 2015) that performs the objective function minimization given in Equation (18.)). [00108] More specific embodiments of the present invention, applied to specific exemplary reaction networks, will now be described with respect to the following two examples; it should be understood that these examples are non-limiting examples: Example 1 - Langmuir Reaction Network

[00109] Assume a reversible reaction between an analyte A and a surface-bound (immobilized) capturing molecule, or ligand, B which is not diffusion or mass transfer limited and obeys pseudo first order kinetics: A + B ¹ AB. (36.)

[00110] The state of this reaction network is given by

{x ₀ (t), x _h (t)} = ([A], [AB])@([B]), (37.) where the concentrations of the analyte A and the product AB are the observed state, and the concentration of the ligand B is the hidden state.

[00111] The differential equations determining the reaction network state over time are

where the analyte concentration [A] is a known constant. [00112] Then the integrated form of Equation (38.b) and Equation

(38.c) are used to eliminate the ligand concentration as

[B] = -[AB] => [B] = R - [AB], (39.) where R is an 'integration constant', namely the initial ligand concentration.

[00113] The elimination steps yield the intermediate differential equation

[AB] - k _a [A](R - [ AB ]) + k _d [AB] = 0, (40.) entirely determined by the observed states but non-linear in the three parameters k _a , k _d and R. This intermediate differential equation defines the discrepancy function.

[00114] The chi-squared objective function of Equation (20.), having said discrepancy function as a parameter, is

where the objective function sums with respect to time the square of the discrepancy function. The chi-squared objective function comprises the discrepancy function defined by Equation (40.) written in vector form.

[00115] The chi-squared expression is reparametrized to the intermediate objective function

com ising a linearized intermediate discrepancy function that is linear with respect to the intermediate parameters k = (k ₁ , k ₂ , k ₃ ). [00116] The intermediate parameters k are obtained with the reparameterization function / _{rep arameterization} defined by

[00117] The reparameterization function has for inverse

[00118] The one or more observed signal(s) (representing observed states within the reaction network), obtained by performing an acquisition method, are used to form an observation matrix as follows:

[00119] The observation matrix for this example, as defined in Equation (24.), is given by wherein [A] is the analyte concentration constant in time and [AB] _t is the product concentration at time t, which is the observed state.

[00120] The observation vector for this example is given by wherein [AB] _t is the derivative against time of the observed signal at time t.

[00121] This example comprises only intermediate parameters shared across all intervals in a sequence of observed signals and no interval specific intermediate parameters. Therefore, all interval observed matrices are defined by Equation (45.) and all interval observed vectors are defined by Equation (46.) for a sequence of observe signals. The observation matrix and observation vector for the sequence of observed signals are given by stacking all interval observation matrices, respectively stacking all interval observation vectors. In the present application "stacking" matrices or vectors, comprises concatenating vertically said matrices or vectors. The intermediate parameters are given by k.

[00122] The first order derivative of the observed signal, which represents the product concentration, is estimated with the central finite difference scheme wherein At is the time interval between two timewise consecutive measurements.

[00123] The values of the intermediate parameters k, which minimize the c ² intermediate objective function, are given by the direct estimator k = -( X ^T X) ^~1 X ^T b , (48.) wherein the intermediate parameters k are entirely determined as an expression of the observation matrix and the observation vector. The diagonal matrix of additional penalty parameters is zero.

[00124] The parameters k _a , k _d and R are then determined by applying the inverse reparameterization function

defined in Equation (44.) to the found values for the intermediate parameters k.

[00125] In this example, we assume Gaussian probability

distributions k^KQ _a R. o and k ₂ ~N(k _d , a ₂ ) (49.) for the intermediate parameters and k ₂ , where the standard deviations and s ₂ have been obtained by an analysis of variance of the values of the discrepancy function. [00126] Then the probability distribution of the estimated R according to the first component in Equation (44.) is the Gaussian ratio distribution which has for expected value

Subtracting the expected value E[R] from twice the estimated value R gives the bias corrected estimate

[00127] In an embodiment, a correction step to obtain the solution to the Langmuir model that minimizes the least square with respect to the observed state can be performed using the method of (Dattner, 2015) that comprises minimizing the objective function defined in Equation (18.).

Example 2 - Langmuir Reaction Network with an Observed Signal Offset

[00128] An observed signal can contain an observed signal offset, which is an offset constant in time by which the observed signal differs from the model signal. Observed signal offsets are particular common in observed signals which have been obtained by performing an acquisition method using a label-free sensor. In this case, to reliably estimate the kinetic parameters of the reaction network, the observed signal offset needs to be jointly determined with the kinetic

parameters. The following embodiment uses a direct estimation method that determines jointly the kinetic parameters and the observed signal offset.

[00129] In such an embodiment the discrepancy function is

[AB] - k _a [A](R - e - [ AB ]) + k _d [AB] + k _d e = 0, (53.) wherein e is an observed signal offset to the observed signal as [AB] ® [AB] + e in Equation (40.), and which is overparametrized by R and e that are exchangeable degrees of freedom.

[00130] The estimation of the parameters k _a , k _d , R and e is

performed under the assumption of minimal absolute observed signal offset e, assumption which resolves the overparameterization and is enforced by appropriately constraining or penalizing e.

[00131] The chi-squared objective function of Equation (20.), with an additional penalty term for the observed signal offset, is given by

where the objective function comprises the discrepancy function in vector form and one additional penalty term comprising one

additional parameter l.

[00132] The intermediate objective function is comprising a linearized intermediate discrepancy function that is linear with respect to the intermediate parameters k = (k ₁ , k ₂ , k ₃ , k ₄ ) and has one additional intermediate penalty term with additional parameter l ₄ .

[00133] The intermediate parameters k and the additional parameter l are obtained with the reparameterization function defined by

[00134] The reparameterization has for inverse

[00135] Typically, the analyte concentration [A] is given a sequence of provided values (i.e. provided by the user)

[Ah ® - ® [A] _t ® - ® [A] _n3 , (58.) each imposed to the reaction network for a time interval T _t . The time interval of the i ^th interval is given by [t _{i :L} , t _{i T} .\. The sequence of analyte concentrations gives rise to a respective sequence of interval observed signals. [00136] The one or more observed signal(s) (representing observed states within the reaction network), obtained by performing an acquisition method, are used to form an observation matrix as follows:

[00137] The interval observation matrix for this example, as defined in Equation (24.), is composed of the following shared interval observation matrix and interval specific observation matrix

wherein [4] _j is the analyte concentration constant in time for the i ^th interval and [AB] _t is the product concentration at time t, which is the observed state.

[00138] The interval observation vector for this example is given by wherein [AB\ _t is the derivative against time of the observed signal at time t. [00139] The intermediate parameters for a sequence of observed signals are given by where the intermediate parameters k _lt k ₂ , k ₃ ) are shared across all observed signal intervals and (k ₄₁ , ... , k _{4 i} , ... , k _{4 ri3} ) are interval specific intermediate parameters.

[00140] The observation matrix and observation vector for a sequence of observed signals are obtained using Equation (26.) from the interval observation matrices and interval observation vectors. [00141] The first order derivative [AB] of the observed signal, which represent the product concentration, is estimated with the central finite difference scheme defined in Equation (47.).

[00142] The values of the intermediate parameters k, which minimize the c ² intermediate objective function, are given by the direct estimator k = -{X ^T X + A) ^~1 X ^T b, (62.) wherein the intermediate parameters k are entirely determined as an expression of the observation matrix, the observation vector, and the diagonal matrix of additional penalty parameters A =

diag(0,0,0, l _{4 1} , ... , l _{4 ί} , ... ,l ₄ . Each interval specific intermediate parameter has a respective additional penalty parameter. The value of the additional intermediate penalty parameters (l _{4 1} , ... ,l _{4 ί} , ... , A _{4 Tl3} ) is pre-defined by the user. In one embodiment, the additional

intermediate penalty parameters are all set to zero. In further embodiments, the user applies methods known in the art to

determine the additional intermediate penalty parameters, methods such as cross-validation or Tikhonov regularization.

[00143] The parameters k _a , k _d , k _t , R and e are then determined by applying the inverse reparameterization function in Equation (57.) to the values of the intermediate parameters k.

[00144] When the intermediate parameters are determined from a sequence of observed signals, each interval specific parameter k _{4 i} determines the interval observed signal offset as e _t = (k ₂ [A\ + /c ₃ ) ^_1 · k _4,i -

[00145] Assuming Gaussian distribution for the intermediate parameters k, the bias correction for R remains unchanged from Equation (52.) in the previous example without observed signal offset. The expected value for the observed signal offset is given by wherein the intermediate parameter k ₄ is as well assumed to have Gaussian distribution {k ₄ ' \k ₄ , a^) around the determined value k ₄ with standard deviation s ₄ determined with an analysis of variance.

The bias corrected observed signal offset e _corrected follows from

Equation (35.).

[00146] In an embodiment, a correction step to obtain the solution to the model that minimizes the least square with respect to the observed state can be performed using the method of (Dattner, 2015) that comprises minimizing the objective function defined in Equation (18.).

Example 3 - Quasi-Steady State Mass Transport Reaction Network

[00147] Common observed signals a represented by reaction networks involving an analyte diffusing to a surface and then reacting with a ligand immobilized at the surface, reaction which is described by the reaction network

A ¹ A _s + B ¹ A _S B. (64.)

Hereafter, [A _s ] is the surface analyte concentration, the concentration of the analyte at the surface where the ligand is immobilized. When the diffusion of the analyte to the surface is slow, the surface analyte concentration becomes time dependent and differs from the fixed analyte concentration [A\, reason the surface analyte concentration needs to be modelled.

[00148] The state of this reaction network is described by

(x ₀ (t), x _h (t)} = ([A], [4 _s B])@([i [ B ]), (65.) where the concentrations of the analyte and the complex are observed states, and the concentrations of the surface analyte and ligand are hidden states.

[00149] The differential equations determining the time evolution of the reaction network state over time are

[A _s ] = [A] - k _t [As ] - k _a [A _s ] [B] + k _d [A _s B] = 0 (a)

[B] = -k _a [A _s ] [B] + k _d [A _s B] 0

[A _S ^' B] = k _a 0y [B] - k _d [A _S B], (c) (66.) where we assume the quasi-steady state approximation [A _s ] = 0 to hold true during the entire duration of the observed signal.

[00150] The Equations (66.b) and (66.c) are then used to eliminate by substitution the ligand concentration as

[B] = -[A _S B] => [B] = R - [A _S B], (67.) where R is an 'integration constant', namely the initial ligand concentration.

[00151] The Equations (66.a) and (66.c) are used to eliminate by substitution the surface analyte concentration as

[A _S ] = [A] - [A _S ^' B]. (68.)

^K t

[00152] The two elimination steps yield the intermediate differential equation

[A _S B] - k _a ([ ] - i A _S ^' B ]) ( R - [. A _S B ]) + k _d [A _s B] = 0, (69.) entirely determined by the concentration of the observed states.

Hereafter, the intermediate differential equation defines the

discrepancy function. [00153] In this example, the chi-squared objective function of Equation (20.) is given by

and comprises as parameter the discrepancy function in vector form and one additional penalty term comprising one additional penalty parameter l.

[00154] The chi-squared intermediate objective function is

comprising the linearized intermediate discrepancy function that is linear with respect to the intermediate parameters k = (k ₁ , k ₂ , k ₃ , k ₄ ) and has one additional intermediate penalty term with additional penalty parameter l ₄ . The intermediate parameter k ₄ needs to be penalized adequately. The optimal value of the addition parameter l ₄ is determined through a known regularization scheme such as

Tikhonov regularization. The value of the additional parameter l ₄ depends on the noise in the observed signal.

[00155] The intermediate parameters k are obtained with the reparameterization function defined by

[00156] The reparameterization has for inverse

[00157] The reparameterization of the additional parameter l in the additional penalty term is given by the identity function

f _r reparameterization a ₄ ) = A. (74.)

[00158] The one or more observed signal(s) (representing observed states within the reaction network), obtained with an acquisition method, are used to form an observation matrix as follows: [00159] The observation matrix as defined in Equation (24.) is given by wherein [A] is the analyte concentration constant in time and [AB\ _t is the product concentration at time t, which represents the observed state.

[00160] The observation vector for this example is given by wherein [A _s B] _t is the derivative against time of the observed signal at time t. [00161] This example comprises only intermediate parameters shared across all intervals in a sequence of observed signals and no interval specific intermediate parameters. Therefore, all interval observed matrices are defined by Equation (75.) and all interval observed vectors are defined by Equation (76.) for a sequence of observe signals. The observation matrix and observation vector for the sequence of observed signals are given by stacking all interval observation matrices, respectively stacking all interval observation vectors. The intermediate parameters are given by k. [00162] Extending the mass transport model with an interval observed signal offset for each interval observed signal in a sequence of observed signals is analogous to Example 2, which extends the Langmuir reaction network with interval observed signal offsets.

[00163] The first order derivative [A _S B] of the product concentration is estimated using the central finite difference scheme in Equation (47.).

[00164] The values of the intermediate parameters k, which minimize the c ² intermediate objective function, are given by the direct estimator k = -(X ^T X + A) ^~1 X ^T b, (77.) wherein the intermediate parameters k are entirely determined as an expression of the observation matrix, the observation vector, and the diagonal penalty matrix A = diag(0,0,0, l ₄ ). The value of the additional penalty parameter l ₄ is herein provided by the experimenter. [00165] The parameters k _a , k _d , k _t and R are then determined from the intermediate parameters k using the inverse reparameterization function in Equation (73.). [00166] The bias introduced by the reparameterization on the parameters k _a , k _d , k _t and R is given by

where the probability distribution p(k ) of the intermediate

parameters k is determined with known methods such as an analysis of variance of the values of the discrepancy function.

[00167] In an embodiment, a correction step to obtain the solution to the model that minimizes the least square with respect to the observed state can be performed using the method of (Dattner, 2015) that comprises minimizing the objective function defined in Equation (18.).

Example 4 - Conformational Change Reaction Network

[00168] Common observed signals are typically represented by reaction networks involving an analyte reacting with a ligand immobilized at the sensor surface to form a product; said product can then undergo a conformational state change, which is a change in the shape of a macromolecule, often induced by environmental factors. The reaction of the analyte with the ligand and the conformational state change reaction are described by the reaction network A + B ¹ AB ¹ AB ^* . (79.)

Hereafter, [A] is the analyte concentration, [B] is the concentration of the ligand immobilized at the sensor surface, [AB] is the concentration of the product from the reaction between the analyte and the ligand, and [ AB ^* ] is the concentration of the conformationally changed state of the product [AB]. The sum [AB] + [AB ^* ] is the observed state, and the concentrations of the ligand B, the product AB, and the product AB*, are hidden states.

[00169] The state of this reaction network is described by

{x ₀ (t), x _h (t)} = ([A], [AB] + [AB ^* ])@( [B], [AB ^* ]), (80.) where the concentrations of the analyte and the sum of the

concentrations of the products are observed states, and the

concentrations of the conformational state (or the product state) and ligand are hidden states.

[00170] The differential equations determining the time evolution of the reaction network state over time are

[B] = -k _a [A] [B] + k _d [AB] (a)

[AB] = k _a [A] [B] - k _d [AB] - k ₂ [AB] + k_ ₂ [AB ^* ] ( b )

[AB ] = k ₂ [AB] - k_ ₂ [AB ^* ]. (c)

(81.) [00171] In a first substitution step, the concentration [B] of the ligand immobilized at the sensor surface is substituted by

[B] = -[AB]-[AB ] and [B] = R - [AB] - [AB ^* ], (82.) in the Equations (81. a, b, c), where R is an 'integration constant', namely the initial ligand concentration.

[00172] In a second substitution step, the hidden product states [AB] and [AB ^* ] in the Equations (81. a, b, c) are substituted by the observed state [C] = [AB] + [AB ^* ], and the combination of hidden states [D] = [AB] - [AB ^* ], using [AB] = ([C] + [D])/2 and [AB ^* ] = ([C] - [D])/2. [00173] Eliminating [D] from the Equations (81. a, b, c), after the first and second substitution steps, yields the intermediate differential equation entirely determined by the concentration of the observed state [C], where the parameter k = k ₂ + k_ ₂ has been introduced for

convenience. Hereafter, the intermediate differential equation defines the discrepancy function.

[00174] The intermediate differential equation in Equation (83.) is a linear second order differential equation, which allows for a vector space ( 4 _{1 i} 4 ₂ ) e IR ² of solutions of the form

where A ₀ e E is a constant determined by the intermediate

differential equation. [00175] In this example, a chi-squared objective function of Equation

(20.) is given by

and comprises as parameter the discrepancy function in vector form.

[00176] The chi-squared intermediate objective function is comprising the linearized intermediate discrepancy function that is linear with respect to the intermediate parameters k =

(kl _> k ₂ , ^3' ^4' ^5)·

[00177] The intermediate parameters k are obtained with the reparameterization function defined by

[00178] The reparameterization has for inverse

[00179] The one or more observed signal(s) (representing observed states within the reaction network), obtained with an acquisition method, are used to form an observation matrix as follows:

[00180] The observation matrix as defined in Equation (24.) is given by

wherein [A] _t is the analyte concentration constant in time for each interval and [C] _t = [ AB] _t + [AB ^* ] _t is the sum of the product

concentration [AB] _t , and the concentration [AB ^* ] _t of the

conformationally changed state of the product, at time t, which represents the observed state.

[00181] The observation vector for this example is given by wherein [C] _t is the derivative against time of the observed signal at time t. [00182] The Equation ... comprises only intermediate parameters shared across all intervals in a sequence of observed signals and no interval specific intermediate parameters. Therefore, all interval observed matrices are defined by Equation (89.) and all interval observed vectors are defined by Equation (90.) for a sequence of observed signals. The observation matrix and observation vector for the sequence of observed signals are given by stacking all interval observation matrices, respectively stacking all interval observation vectors. The intermediate parameters are given by k.

[00183] The derivatives [c] and [C] of the observed states are estimated using an interpolation of the data as given in Equation (29.).

[00184] The values of the intermediate parameters k, which minimize the c ² intermediate objective function, are given by a direct estimator k = -(X ^T X + A)- ¹ X ^T b (91.) wherein the intermediate parameters k are entirely determined as an expression of the observation matrix, and the observation vector. [00185] Due to the vector space of solutions to the intermediate differential equation of the conformational change model, vector space of solutions given by Equation (84.), the direct estimator in Equation (91.) is degenerate, yielding a vector space of solutions for k with equal value of the objective function . Within the vector space of solutions for k, there exists a unique solution that minimizes as well the additional constraints {f tiai {P’ ^e i)’- , fS _itial {p, <4*)} for the initial conditions of the intermediate differential equation.

[00186] A first initial condition for the differential equation is obtained from the differential Equation (81. b) at time zero, given by

[00187] A second initial condition for the differential equation is obtained the differential Equations (81.) at equilibrium, and is given This initial condition can be used when one of the intervals reached the equilibrium state.

[00188] The parameters k _a , k _d , k ₂ , k_ ₂ and R are then determined from the intermediate parameters k by the minimization

which can be performed using standard optimization methods such as Levenberg-Marquardt. This optimization has a very low computational cost compared to fitting the complete model to all data points. [00189] In an embodiment, a correction step to obtain the solution to the model that minimizes the least square with respect to the observed state can be performed using the method of (Dattner, 2015) that comprises minimizing the objective function defined in Equation (18.)·

[00190] In the above-mentioned examples there are described steps in which the intermediate objective function is determined. However, in a further embodiment of the present invention the intermediate objective function is predetermined, so said above-described steps of determining the intermediate objective function are not performed. For example, there may be provided a library or memory, containing a plurality of intermediate objective functions, each for a different reaction network. In such a case the present invention may involve selecting an intermediate objective function from the library or memory.

[00191] Additionally, in a further aspect of the present invention there is provided (preferably on a tangible data carrier) a program code arranged for causing a processor to carry out the steps described above, when said processor executes said program code to determine the kinetic parameters of the reaction network.

[00192] In one embodiment the program code may be configured to receive a user's selection of an intermediate objective function from the library or memory; and then to carry out the steps described above using said selected intermediate objective function to determine the kinetic parameters of the reaction network.

[00193] Various modifications and variations to the described embodiments of the invention will be apparent to those skilled in the art without departing from the scope of the invention as defined in the appended claims. Although the invention has been described in connection with specific preferred embodiments, it should be understood that the invention as claimed should not be unduly limited to such specific embodiment.