You have a mess of parenthesis. I think your system should look like:
Clear[t];
a = 0.00001381594775713100;
c = 0.00007261477318255780;
j = 0.00004126704565565830;
d = 0.00006126704565565830;
b = 0.498812483;
r = 0.498812483;
Population = 40000;
mi = 0.4;
Na = Population*mi;
Ns = Population*(1 - mi);
ts = 400;
system = {
Iz'[5t]/5 == a Iz[5t](Na - Iz[5t]) + j Iz[5t](Ns - Is[5t]) - b Iz[5t],
Is'[t] == c Is[t](Ns - Is[t]) + d Is[t](Na - Iz[t]) - r Is[t]
};
initialvalues = {Is[0] == 160, Iz[0] == 40};
solution = NDSolve[Join[system, initialvalues], {Iz, Is}, {t, 0, ts}];
But then it is a sort of Delay Differential Equation. So if you run it you get an error:
Initial history needs to be specified for all variables for delay-differential equations. >>
And if you look at the help DDLs have usually this sort of initial conditions:
ndssol = First[
NDSolve[{x'[t] == y[t], y'[t] == -x[t - 1], x[t /; t <= 0] == t^2,
y[t /; t <= 0] == 2 t}, {x, y}, {t, 0, 5}]];
Plot[Evaluate[{x[t], y[t]} /. ndssol], {t, 0, 5}]
But then if here for example you change time shift into time rescaling like:
ndssol = First[
NDSolve[{x'[t] == y[t], y'[t] == -x[2 t], x[t /; t <= 0] == t^2,
y[t /; t <= 0] == 2 t}, {x, y}, {t, 0, 5}]];
Plot[Evaluate[{x[t], y[t]} /. ndssol], {t, 0, 5}]
You get an error:
The method currently implemented for delay differential equations does not support delays that depend directly on the time variable or dependent variables. >>
So I think we might be in trouble here. An expert opinion is needed.