  Galaxy evolution emerges from the balance between cosmic gas accretion,
fueling star formation, and supernova (SN) feedback, regulating the metal
enrichment. Hence, the stellar mass ($M_*$) - gas metallicity relation (MZR) is
key to understand the physics of galaxies. High-quality JWST data enable
accurate measurements of the MZR up to redshift z=10. Our aims are to
understand the observed MZR, its connection with the star formation rate (SFR),
the role played by SFR stochasticity, and how it is regulated by SN feedback.
We compare the MZR from the JADES, CEERS, and UNCOVER surveys, which comprise
about 180 galaxies at $z=3-10$ with $10^6<M_*/M_\odot<10^{10}$, with 200
galaxies from the SERRA cosmological simulations. To interpret the MZR, we
develop a minimal model for galaxy evolution that includes: cosmic accretion
modulated with an amplitude $A_{100}$ on 100 Myr; a time delay $t_d$ between
SFR and SN; SN-driven outflows with a varying mass loading factor
$\epsilon_{SN}$. Using our minimal model, we find the observed mean MZR is
reproduced by weak outflows ($\epsilon_{SN}=1/4$), in line with findings from
JADES. Matching the observed MZR dispersion requires $t_d=20$ Myr and a
$A_{100}=1/3$ modulation of the accretion rate. Successful models have low
stochasticity ($\sigma_{SFR}=0.2$), yielding a MZR dispersion of
$\sigma_{Z}=0.2$. Such values are close but lower than SERRA predictions
($\sigma_{SFR}=0.24$, $\sigma_{Z}=0.3$), clarifying why SERRA show no clear MZR
trend and some tension with the observations. As the MZR is very sensitive to
SFR stochasticity, models predicting high r.m.s. values ($\sigma_{SFR}=0.5$)
result in a ``chemical chaos'' (i.e. $\sigma_{Z}=1.4$), virtually destroying
the MZR. As a consequence, invoking a highly stochastic SFR
($\sigma_{SFR}=0.8$) to explain the overabundance of bright, super-early
galaxies leads to inconsistencies with the observed MZR.
