data = {{18.283, 0.0003365}, {39.415, 0.0003892}, {60.547,
    0.00045}, {81.679, 0.0005631}, {102.811, 0.0006446}, {123.943,
    0.0006756}, {145.075, 0.0007655}, {166.207, 0.0008306}, {187.339,
    0.0008814}, {208.471, 0.000896}, {229.603, 0.0009816}, {250.735,
    0.0011191}, {271.867, 0.0010797}, {292.999, 0.0011473}, {314.131,
    0.0012265}, {335.263, 0.0013187}, {356.395, 0.0012809}, {377.527,
    0.0013523}, {398.659, 0.0013841}, {419.791, 0.0014507}, {440.923,
    0.0014885}, {462.055, 0.0014898}, {483.187, 0.0015491}, {504.319,
    0.0015838}, {525.451, 0.0015651}, {546.583, 0.0015903}, {567.715,
    0.0016328}, {588.847, 0.0016423}, {609.979, 0.0016149}, {631.111,
    0.0016508}, {652.243, 0.0016296}, {673.375, 0.0016465}, {694.507,
    0.0016776}, {715.639, 0.0016668}, {736.771, 0.0016685}, {757.903,
    0.0016967}, {779.035, 0.001717}, {800.167, 0.0016997}, {821.299,
    0.001721}};
tmax = Max[data[[All, 1]]];
e0 == 0.0002;
k21 == 0.001;
kh == 0.0001;
model = ParametricNDSolveValue[{s'[
      t] == -16  s[t] (0.0002 - x1[t] - x2[t]) + k21 x1[t] - kh  s[t],
    x1'[t] == 16  s[t] (0.0002 - x1[t] - x2[t]) - (k21 + k2) x1[t],
    x2'[t] == k2 x1[t] - k3 x2[t],
    p1'[t] == k2 x1[t] + kh s[t],
    s[0] == 0.002, x1[0] == 0, x2[0] == 0, p1[0] == 0},
   p1, {t, 0, tmax}, {k2, k3}];
fit = NonlinearModelFit[data,
  model[k2, k3][t], {{k2, 0.4}, {k3, 0.7}}, {t, 0, tmax}]
plotfit = Plot[fit[t], {t, 0, tmax}, PlotRange -> {0, 0.0021}];
plotdata = ListPlot[data, PlotStyle -> PointSize[0.01]];
Show[plotfit, plotdata]

NDSolve::deqn: Equation or list of equations expected instead of True in the first argument
General::ivar: 0 is not a valid variable. >>
ParametricNDSolve::dsvar: 0.016777965285714284` cannot be used as a variable. >>
General::ivar: 0.016777965285714284` is not a valid variable. >>