I have written a code in Mathematica for a simultaneous system of ordinary differential equations. Though the fitting is just ok, but I want to make it more accurate. I am attaching my code here. Could you please give me a suggestion on how to make my model fitting better? Please note I am a new user of Mathematica, and I have prepared this code already with help of a user of this community, so if you suggest a change please clearly mention/highlight it so that I can easily understand it. I want to use similar fitting procedure as already in code, do I need to add another ODE equation or what? I am interested in knowing how can I describe ODE for 'SD' such that it flatten out not at zero but somewhere close to real data {60, 0.003234}.

bcdata = {{0, 0.000836}, {1, 0.001650}, {2, 0.002068}, {3, 0.002772}, {4, 0.004730}, {5, 0.005390}, {10, 0.006248}, {20, 0.006336}, {40, 0.006512}, {60, 0.006534}};
aqdata = {{0, 0.003564}, {1, 0.003630}, {2, 0.004818}, {3, 0.005786}, {4, 0.007260}, {5, 0.007964}, {10, 0.008998}, {20, 0.008646}, {40, 0.008492}, {60, 0.008316}};
gasdata = {{0, 0.000000}, {1, 0.000352}, {2, 0.000638}, {3, 0.001188}, {4, 0.001474}, {5, 0.001650}, {10, 0.002684}, {20, 0.003564}, {40, 0.003564}, {60, 0.003916}};
sddata = {{0, 0.017622}, {1, 0.016368}, {2, 0.014476}, {3, 0.012254}, {4, 0.008536}, {5, 0.006996}, {10, 0.004048}, {20, 0.003454}, {40, 0.003432}, {60, 0.003234}};
eqns = {bc'[t] == k1 sd[t] - k4 bc[t], aq'[t] == k2 sd[t] - k5 aq[t], gas'[t] == k3 sd[t] + k4 bc[t] + k5 aq[t], sd'[t] == (-k1 - k2 - k3) sd[t]};
Thread[{aa[t_], bb[t_], cc[t_], dd[t_]} = DSolveValue[{eqns, bc[0] == 0.000836, aq[0] == 0.003564, gas[0] == 0.000000, sd[0] == 0.017622}, {bc[t], aq[t], gas[t], sd[t]}, t]];
model[k1_, k2_, k3_, k4_, k5_] := Sum[(aa[bcdata[[i, 1]]] - bcdata[[i, 2]])^2 + (bb[aqdata[[i, 1]]] - aqdata[[i, 2]])^2 + (cc[gasdata[[i, 1]]] - gasdata[[i, 2]])^2 + (dd[sddata[[i, 1]]] - sddata[[i, 2]])^2, {i, Length@bcdata}]
fit = Last@NMinimize[model[k1, k2, k3, k4, k5], {k1, k2, k3, k4, k5}]
Thread[{k1, k2, k3, k4, k5} = Values@fit];
Show[Plot[{aa[t], bb[t], cc[t], dd[t]}, {t, 0, 60}, Frame -> True], ListPlot[{bcdata, aqdata, gasdata, sddata}]]