Abstract Established genome-scale modelling methods primarily predict reaction fluxes while established high-throughput experimental technologies primarily measure molecular species concentrations. This paradoxical situation has arisen because of the problem to implement the nonlinear constraints that represent reaction kinetic rate equations without recourse to convenient yet inaccurate approximations. We present a mathematically elegant and computationally tractable solution to this problem. First we introduce a mathematical reformulation of established knowledge on metabolic reactions and reaction kinetics in matrix vector notation. Then we introduce variational kinetics, a novel approach to satisfy steady-state reaction kinetics at genome scale using novel mathematical and numerical optimisation techniques. We illustrate how this approach may be used to simultaneously optimise over the set of steady-state reaction fluxes, thermodynamically feasible kinetic parameters and kinetically feasible elementary and phenomenological rate laws, with solution times competitive with linear optimisation.