Improving bad fit for data using Matlab

Question

I am trying to fit some data to estimate model parameters in a system of 5 ODEs. I cannot get a good fit. First I tried using Matlab lsqcurvefit and I get a very poor fit.
Then I tried fmincon and this is the code that I used,

    close all
clear
initial= [1.8;0.55;8.5*10^-12;0.095;0.25;10^6;0.35];

t=[0.509283106;1.564581803;2.607620249;3.6594436;4.706343542;5.708739739;6.695689949;8.810390184;9.867523091;...
    10.86417531;11.92179091;12.93253759;13.97605872;];

data=[3.288981289;4.875259875;4.613305613;5.675675676;5.995841996;6.607068607;4.88981289;5.180873181;0.043659044...
    ;3.288981289;5.224532225;3.594594595;3.405405405;];

lb=[1,0.2,10^-13,0.001,0.001,10^5,0.001];
ub=[100,100,inf,20,20,10^9,20];

x=[t,data];
[fit,fval]=fmincon(@(parameters)ssq(parameters,x),initial,[],[],[],[],lb,ub);
figure
scatter(t,data,5,'r')
hold on

time=0:14;
fittedValues=model(fit',time);
plot(time,fittedValues,'g')
legend('data','fit')

function SSE=ssq(parameters,x)
time=x(:,1);
data=x(:,2);

fitted=model(parameters,time);
error=fitted(:,1)-data;
SSE=sum(error.^2);

end

function output= model(parameters,thours)
    x0=[10^3,0.001,0.001,10^-8];
    [time,values] = ode45(@Equations,thours,x0);

function s=Equations(t,y)

a1=0.3368;
p8=0.03;
p2=10^-3;
N=(((2.5+7.5)/2)*10^9)*5;
p8=1/24;
a2=1/11;
k2=10^6/a2;
p5=0.75;
p6=8;
p7=0.15;

s=zeros(4,1);

s(1)=(parameters(4)*y(1))*(1-(y(1)/parameters(6)))-p7*a1*y(1)*(1-(y(3)/k2))-...
    parameters(1)*y(4)*y(1)/(parameters(2)*y(4)+y(1))-p2*y(1)...
    +p5*y(3)+p8*y(2);
s(2)=p8*parameters(1)*y(4)*y(1)/(parameters(2)*y(4)+y(1))-p8*y(2)+...
    parameters(7)*y(2)*(1-(y(2)/(y(4)*p6)));
s(3)=p7*a1*y(1)*(1-(y(3)/k2))-p5*y(3)+parameters(5)*y(3);
s(4)=parameters(3)*(N-y(4))*y(1)-p8*y(4);

end


output=values(:,1)+values(:,2)+values(:,3);
output=log10(output);

end

With this code I had to increase the upper bound of 3rd parameter to inf as I kept on geting warnings as, Derivative finite-differencing step was artificially reduced to be within bound constraints. This may adversely affect convergence. Increasing distance between bound constraints, in dimension 3, to be at least 1.9983e-08 may improve results.
The fmincon fit is

I believe the fits depends highly on what I choose as the initial guesses and the lower and upper bounds. What can I do to get a much better fit?
Because when I changed the model parameters by hand (just arbitrarily changing parameter values and plotting the model output against data) I seem to get get somewhat of a better ft than I am getting using these data fitting options.

So, can someone please point to me what I am doing wrong here and some ideas as to how I can achieve a better fit.