Coupled parabolic differential equations with time delay

Is it possible for NDSolve to solve delay partial differential equations with simple Neumann boundary conditions?

An example I have is as below:

ClearAll["Global*"]
A1 = 350; A2 = 100;
μ1 = 1/100; α1 = 5/100; β1 =
5*10^-6; a2 = 100; β2 = 7*10^-7; μ2 =
391*10^-9; α2 = 3/10; γ = 1/100; D1 = 2; D2 = 1;
τ1 = 1/5;
τ2 = 1/6;
omega = ImplicitRegion[0 <= x <= 3, x];
sol = NDSolve[{D[s1[t, x], t] - D1*Laplacian[s1[t, x], x] -
A1 - μ1*s1[t, x] - β1*s1[t, x]*I1[t, x] ==
NeumannValue[0, True],
D[I1[t, x], t] -
D1*Laplacian[I1[t, x], x] - β1*Exp[-μ1*τ1]*
s1[t - τ1, x]*I1[t - τ1, x] + (μ1 + α1)*
I1[t, x] == NeumannValue[0, True],
D[s2[t, x], t] - D2*Laplacian[s2[t, x], x] -
A2 - μ2*s2[t, x] - β2*s2[t, x]*I2[t, x] ==
NeumannValue[0, True],
D[I2[t, x], t] -
D2*Laplacian[I2[t, x], x] - β2*Exp[-μ2*τ2]*
s2[t - τ2, x]*
I1[t - τ2, x] + (μ2 + α2 + γ)*I2[t, x] ==
NeumannValue[0, True],
D[R2[t, x], t] -
D2*Laplacian[I2[t, x], x] - γ*I2[t, x] + μ2*R2[t, x] ==
NeumannValue[0, True],
s1[0, x] == 11 + Sin[2 x],
I1[t /; t <= 0, x] == 1/10 + 1/100 Sin[2 x],
s2[0, x] == 100 + Sin[2 x], I2[t /; t <= 0, x] == 0, R2[0, x] == 0
}, {s1[t, x], I1[t, x], s2[t, x], I2[t, x], R2[t, x]}, {t, 0, 300},
Element[x, omega]]


But there is no useful outputs.

From the documents of NDSolve, there are only illustrative examples for delay ordinary differential equations.

• The numerical solutions to all the dependent variables, s1,s2,I1,I2 and R2, in interpolated function forms as NDSolve usually does, are expected outputs. Since this is a linear system, exact solutions may also be possible. – LCFactorization May 6 '17 at 3:26
• Laplacian operator is readily for higher dimensional extensions – LCFactorization May 6 '17 at 3:29
• NDSolve::delpde: Delay partial differential equations are not currently supported by NDSolve. – zhk May 6 '17 at 3:31
• thank you. I am not sure whether such feature will be available in the future releases. Delay pdes are still hot topics in recent publications. So it seems in order to use NDSolve` to solve the system, one has to do a lot of basic job manually – LCFactorization May 6 '17 at 3:36
• Maybe this is of some kind of help, scicomp.stackexchange.com/questions/23402/… – zhk May 6 '17 at 3:44