kp1 = 1;
km1= 0.1;
k2 = 0.01;

sol = NDSolve[{
e'[t] == -kp1 e[t] s[t] + km1 es[t] + k2 es[t], 
s'[t] == -kp1 e[t] s[t] + km1 es[t], 
es'[t] == kp1 e[t] s[t] - (km1 + k2) es[t], 
p'[t] == k2 es[t], 
e[0] == 0.01, 
s[0] == 1.0, 
es[0] ==  0, 
p[0] == 0}, {e, s, p, es}, {t,  50000}]

Plot[Evaluate[{s[t]/.sol,p[t]/.sol, 10 e[t]/.sol, 10 es[t]/.sol}], {t, 0, 50000}, PlotRange -> All,AxesLabel->{"t", "Concentration"},  PlotLegends->{"s(t)", "p(t)", "e(t) × 10", "es(t) × 10"}]

g1 = ParametricPlot[Evaluate[{s[t], p'[t]}/.sol], {t, 0, 50000}, AspectRatio -> 1, PlotRange -> {{0,1},{0, 10^(-4)}},AxesLabel->{"s(t)", "p'(t)"}]

vmax = k2 e[0]/.sol

km = (km1 + k2)/kp1

v0[substrat_] := vmax substrat/ (km + substrat)

g2 = Plot[v0[substrat], {substrat, 0, 1}, PlotRange -> {{0, 1},{0, 10^(-4)}}, AspectRatio -> 1, AxesLabel->{"s(t)", "p'(t) ≡ (Subscript[v, 0])"}, PlotStyle -> {Orange,Dashed}]

Show[g1, g2]