Error when using NDSolveValue. Lists are not the same shape

Please, help me with the following issue.
I am solving the following system of the coupled PDEs

with the following boundary conditions

and initial conditions.

After performing the following simplified code (for simplicity I set all parameters equal to 1):

(*specify the system*)
op = 
  {-Inactive[Grad][n[t, x], t] + Laplacian[n[t, x], x] - 
     Inactive[Grad][n[t, x]*Grad[a[t, x], x], x] + a[t, x]*b[t, x]  -  
     n[t, x] - (n[t, x])^2,
   -Inactive[Grad][a[t, x], t] + Inactive[Laplacian][a[t, x],x] + 
     (1 + Tanh[1 - b[t, x]] ) - (1 + b[t, x])*a[t, x],
   -Inactive[Grad][b[t, x], t] + 
      Inactive[Grad][n[t, x]*Inactive[Grad][b[t, x], x], x] - 
      Inactive[Grad][n[t, x], x] + n[t, x]*Inactive[Grad][a[t, x], x]}

(*specify the Dirichlet boundary conditions*)
bcs = 
  {DirichletCondition[n[t, x] == E^-t, x == 0], 
   DirichletCondition[b[t, x] == 1, x == 0]};

(*specify the Neumann boundary conditions (for function a(t,x))*)
BNeuma = {NeumannValue[-a[t, x], x == 0], NeumannValue[0, x == 1]};
BNeumn = NeumannValue[0, x == 1];
BNeumb = NeumannValue[0, x == 1];
 
(*function to use in the picewise intial conditions*)
xTild = 0.05;
n0[x_] := 1/xTild^3*(x - xTild)*(2*x^2 - xTild*x - xTild^2);
b0[x_] := 1/xTild^3 (x - xTild)*(2*x^2 - xTild*x - xTild^2);

(*initial conditions themselfs*)
incs = 
 {n[0, x] == Piecewise[{{n0[x], x < xTild}, {0, x > xTild}}], 
  a[0, x] == 0, 
  b[0, x] == Piecewise[{{b0[x], x < xTild}, {0, x > xTild}}]};

(*solve the system*)
{nfun, afun, bfun} = 
  NDSolveValue[
    {op == {BNeumn, BNeuma, BNeumn}, bcs, incs}, 
    {n, a, b}, {x, 0, 1}, {t, 0, 10}];

I get the following error:

Set::shape: Lists {nfun, afun, bfun} and NDSolveValue[<<1>>] are not the same shape.

Please, give me a hint, why my code is incorrect. If it is nessesary, I can upload a full code with all parameters, but it gives exactly the same error. My Mathematica version is 12.1.

eq = {-D[n[t, x], t] + D[n[t, x], x, x] - 
     D[n[t, x]*D[a[t, x], x], x] + a[t, x]*b[t, x] - 
     n[t, x] - (n[t, x])^2 == 
    0, -D[a[t, x], t] + 
     D[a[t, x], x, x] + (1 + Tanh[1 - b[t, x]]) - (1 + b[t, x])*
      a[t, x] == 
    0, -D[b[t, x], t] + D[n[t, x]*D[b[t, x], x], x] - D[n[t, x], x] + 
     n[t, x]*D[a[t, x], x] == 0};

(*specify the Dirichlet boundary conditions*)
bcs = {n[t, x] == E^-t /. x -> 0, 
   b[t, x] == 1 /. 
    x -> 0, -Derivative[0, 1][a][t, x] + a[t, x] == 0 /. x -> 0, 
   Derivative[0, 1][n][t, 1] == 0, Derivative[0, 1][a][t, 1] == 0, 
   Derivative[0, 1][b][t, 1] == 0};

(*specify the Neumann boundary conditions (for function a(t,x))*)
BNeuma = NeumannValue[-a[t, x], x == 0];

(*function to use in the picewise intial conditions*)
xTild = 0.05;
n0[x_] := 1/xTild^3*(x - xTild)*(2*x^2 - xTild*x - xTild^2);
b0[x_] := 1/xTild^3 (x - xTild)*(2*x^2 - xTild*x - xTild^2);

(*initial conditions themselfs*)
incs = {n[0, x] == If[x <= xTild, n0[x], 0], a[0, x] == 0, 
   b[0, x] == If[x < xTild, b0[x], 0]};


{nfun, afun, bfun} = 
  NDSolveValue[Join[eq, bcs, incs], {n, a, b}, {x, 0, 1}, {t, 0, 10}];

{DensityPlot[nfun[t, x], {x, 0, 1}, {t, 0, 10}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> Automatic, PlotRange -> All], 
 DensityPlot[afun[t, x], {x, 0, 1}, {t, 0, 10}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> Automatic, PlotRange -> All], 
 DensityPlot[bfun[t, x], {x, 0, 1}, {t, 0, 10}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> Automatic, PlotRange -> All]}