Neuman boundary conditions for system of PDE in 2D

The system of equations for three variables: n(x,y,t), b(x,y,t) and a(x,y,t) System

The rest of letters are coefficients. System is defined over a square region with the lower left corner at the origin.

Boundary conditions: Boundary

Bold n is the normal vector to the boundary (Neuman condition).

And initial conditions: Initial

My code for solution is the following

(*region*) \[CapitalOmega] = Rectangle[{0, 0}, {1, 1}];  (*constants*) \[Mu]n = 0.001; \[Mu]b = 0.001; \[Chi] = 0.1; bHat = 1.5; \[Lambda]0 = 100; \[Lambda]1 = 100; \[Lambda]2 = 10; \[Lambda]4 = 100; \[Lambda]5 = 10; nHat = 1; \[Lambda]7 = 10; \[Delta] = 0.01; \[Alpha] = 2.5;  aNeuman =    NeumannValue[\[Lambda]7*bHat*a[t, x, y], {x, y} \[Element]RegionBoundary[\[CapitalOmega]]];  (*equations*) eq = {-D[n[t, x, y], t] + \[Mu]n*Laplacian[n[t, x, y], {x, y}] - \[Chi] *       Div[n[t, x, y]*Grad[a[t, x, y], {x, y}], {x, y}] + \[Lambda]1*       a[t, x, y]*b[t, x, y] - \[Lambda]2*n[t, x, y] - \[Lambda]0*(n[t, x, y])^2 == 0,       -D[a[t, x, y], t] + Laplacian[a[t, x, y], {x, y}] +        \[Lambda]4 /2*(1 + Tanh[(1 - b[t, x, y])/\[Delta]]) - (\[Lambda]4 + b[t, x, y])*a[t, x, y] == aNeuman,        -D[b[t, x, y], t] + \[Mu]b*Div[n[t, x, y]*Grad[b[t, x, y], {x, y}], {x, y}] -       Norm[\[Mu]n*Grad[n[t, x, y], {x, y}] - \[Chi] *n[t, x, y]*Grad[a[t, x, y], {x, y}]] == 0};  (*initial conditions*) incs = {n[0, x, y] == If[x == 0 || y == 0 || x == 1 || y == 1, nHat, 0],         a[0, x, y] == 0,         b[0, x, y] == If[x == 0 || y == 0 || x == 1 || y == 1, bHat, 0]};  (*boundary conditions*) bcs = {DirichletCondition[n[t, x, y] == nHat*E^(-\[Alpha]*t), True],        DirichletCondition[b[t, x, y] == bHat, True]};  (*solution*) {nfun, afun, bfun} = NDSolveValue[Join[eq, bcs, incs], {n, a, b}, {t, 0, 2}, {x, y} \[Element] \[CapitalOmega]]; 

This code gives the error (apparently because of the wrong way of adding Neumann condition)

LinearSolve::parpiv: Zero pivot was detected during the numerical factorization  or there was a problem in the iterative refinement process.  It is possible that the matrix is ill-conditioned or singular. 

I actually have two questions.

  1. How to properly implement the boundary condition for a(x,y,t)?

  2. When I replace this Neuman condition with the Dirichlet one (which is wrong, but for testing puprposes), the solver starts to work, but then I get another error:

    NDSolveValue::ndsz: At t == 1.8768713852597982`, step size is effectively zero; singularity or stiff system suspected.

This error goes away when I set all coeffitienst of the equations equal to 1. But could you, please, help me with solving the system with original coeffitients?

P.S. My Mathematica version is 12.1

P.P.S. Sorry for that Greek letters. I tried make them consistent with the original equations. If they really distract you, please, tell, I will rename.