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.

Sam,

I´m not sure if this system has a kind of delay, despite the error message from Mathematica.

Actually, when we think about a system with delay, all equations have the same time scale for the derivatives. In others words, the left side of equations are in function of "t", and in the rigth side we have something in "t - x" or "t + x".

Thanks.