We investigate correlations between different physical properties of star-forming galaxies in the Evolution and Assembly of GaLaxies and their Environments (EAGLE) cosmological hydrodynamical simulation suite over the redshift range $0le zle 4.5$. A principal component analysis reveals that neutral gas fraction ($f_{rm gas, neutral}$), stellar mass ($M_{rm stellar}$) and star formation rate (SFR) account for most of the variance seen in the population, with galaxies tracing a two-dimensional, nearly flat, surface in the three-dimensional space of $f_{rm gas, neutral}-M_{rm stellar}-rm SFR$ with little scatter. The location of this plane varies little with redshift, whereas galaxies themselves move along the plane as their $f_{rm gas, neutral}$ and SFR drop with redshift. The positions of galaxies along the plane are highly correlated with gas metallicity. The metallicity can therefore be robustly predicted from $f_{rm gas, neutral}$, or from the $M_{rm stellar}$ and SFR. We argue that the appearance of this fundamental plane of star formation is a consequence of self-regulation, with the planes curvature set by the dependence of the SFR on gas density and metallicity. We analyse a large compilation of observations spanning the redshift range $0lesssim rm zlesssim 2.5$, and find that such a plane is also present in the data. The properties of the observed fundamental plane of star formation are in good agreement with EAGLEs predictions.