Need help!!! I am modeling the reaction kinetics of enzyme mechanism of a reaction system. Based on the reaction rate equations, 8 differential equations with 6 rate constants (k1-k6) are generated. The results of the expected response variable "b" are available as provided in the code below. I want to fit the model to the experimental data (b), in order to come up with the suitable rate constants for the reaction system. I have tried to adopt my problem to the solution provided by IIian but could not get reasonable results for k1-k6 and meaningful chart. Your help will be highly appreciated. Thank you.

Clear[k1, k2, k3, k4, k5, k6]; 

totaltime = 6; 

dataset1 = {{0, 0}, {5, 0.5472}, {10, 0.6310}, {15, 0.7120}, {20, 0.7630}, {25, 0.846992}, {30, 0.8364}}; 

model = ParametricNDSolveValue[{s'[t] == -k1*s[t]*e[t] + k2*x[t], 

a'[t] == -k4*a[t]*e1[t] + k5*y[t], 

b'[t] == k3*x[t], 

g'[t] == k6*y[t], 

e'[t] == -k1*s[t]*e[t] + k2*x[t] + k6*y[t], 

e1'[t] == -k4*s[t]*e[t] + k5*y[t] + k3*x[t], 

x'[t] == k1*s[t]*e[t] - k2*x[t] - k3*x[t], 

y'[t] == k4*a[t]*e1[t] - k5*y[t] - k6*y[t], 

s[0] == 4, a[0] == 0.91, e[0] == 0.36, b[0] == 0, g[0] == 0, e1[0] == 0, x[0] == 0, y[0] == 0}, y, {t, 0, totaltime}, {k1, k2, k3, k4, k5, k6}]; 

fit = FindFit[dataset1, model[k1, k2, k3, k4, k5, k6][t], {k1, k2, k3, k4, k5, {k6, 0}}, t] 

Plot[model[k1, k2, k3, k4, k5, k6][t] /. fit, {t, 0, 30}, Epilog -> {Red, Point[dataset1]}]