  Chemical reaction networks offer a natural nonlinear generalisation of linear
Markov jump processes on a finite state-space. In this paper, we analyse the
dynamical large deviations of such models, starting from their microscopic
version, the chemical master equation. By taking a large-volume limit, we show
that those systems can be described by a path integral formalism over a
Lagrangian functional of concentrations and chemical fluxes. This Lagrangian is
dual to a Hamiltonian, whose trajectories correspond to the most likely
evolution of the system given its boundary conditions. The same can be done for
a system biased on time-averaged concentrations and currents, yielding a biased
Hamiltonian whose trajectories are optimal paths conditioned on those
observables. The appropriate boundary conditions turn out to be mixed, so that,
in the long time limit, those trajectories converge to well-defined attractors.
We are then able to identify the largest value that the Hamiltonian takes over
those attractors with the scaled cumulant generating function of our
observables, providing a non-linear equivalent to the well-known
Donsker-Varadhan formula for jump processes. On that basis, we prove that
chemical reaction networks that are deterministically multistable generically
undergo first-order dynamical phase transitions in the vicinity of zero bias.
We illustrate that fact through a simple bistable model called the Schl\"ogl
model, as well as multistable and unstable generalisations of it, and we make a
few surprising observations regarding the stability of deterministic fixed
points, and the breaking of ergodicity in the large-volume limit.
