  We apply a Gaussian variational approximation to model reduction in large
biochemical networks of unary and binary reactions. We focus on a small subset
of variables (subnetwork) of interest, e.g. because they are accessible
experimentally, embedded in a larger network (bulk). The key goal is to write
dynamical equations reduced to the subnetwork but still retaining the effects
of the bulk. As a result, the subnetwork-reduced dynamics contains a memory
term and an extrinsic noise term with non-trivial temporal correlations. We
first derive expressions for this memory and noise in the linearized (Gaussian)
dynamics and then use a perturbative power expansion to obtain first order
nonlinear corrections. For the case of vanishing intrinsic noise, our
description is explicitly shown to be equivalent to projection methods up to
quadratic terms, but it is applicable also in the presence of stochastic
fluctuations in the original dynamics. An example from the Epidermal Growth
Factor Receptor (EGFR) signalling pathway is provided to probe the increased
prediction accuracy and computational efficiency of our method.
