I am trying to solve the pde (Eqn (4)) in the following paper with the functions given in the methods section. I want to have plots as in Fig 2b,c.

https://home.gwu.edu/~glan/files/2012.Nature.Physics%20-%20The%20energy%A8Cspeed%A8Caccuracy%20trade-off%20in%20sensory%20adaptation.All.pdf

I am using NDSolve for it, but I get an error that says:

At t == 0.0021905285143251785`, step size is effectively zero; singularity or stiff system suspected.

The paper doesn't have any initial condition so I constructed one that satisfies the boundary conditions (based on a stochastic simulation). Here is my code:

ClearAll["Global'*"]
K0 = 1;
Kd[m_] = K0 E^(2 m);
H = 1;
s = 20;
G[m_] = 1/(1 + (s/Kd[m])^H);
\[Omega]m = 5;
\[Omega]a = 50;
\[CapitalDelta]m =  (0.01)*5;
\[CapitalDelta]a = (0.01)*\[Omega]a;
Fa[a_, m_] = -\[Omega]a*(a - G[m]);
a0 = 0.5;
C1 = (\[CapitalDelta]m \[Omega]a)/(\[CapitalDelta]a \[Omega]m);
\[Beta] = 1;
Fm[a_, m_] = -\[Omega]m*(a - a0)*(\[Beta] - (1 - \[Beta])*C1*\!\(
\*SubscriptBox[\(\[PartialD]\), \(m\)]\(G[m]\)\));
f[a_, m_] =  \!\(
\*SubscriptBox[\(\[PartialD]\), \(a\)]\(Fa[a, m]\)\) - \!\(
\*SubscriptBox[\(\[PartialD]\), \(m\)]\(Fm[a, m]\)\);
F[a_, m_] = {Fa[a, m], Fm[a, m]};
M = {{\[CapitalDelta]a, 0}, {0, \[CapitalDelta]m}};
eqn = D[P[a, m, t], t] + f[a, m]*P[a, m, t] + 
   Dot[F[a, m], Grad[P[a, m, t], {a, m}]] + 
   Div[Dot[-M, Grad[P[a, m, t], {a, m}]], {a, m}] ;
T = 10;
nv1 = NeumannValue[-Fm[a, m]/\[CapitalDelta]m P[a, m, t], m == 0];
nv2 = NeumannValue[-Fm[a, m]/\[CapitalDelta]m P[a, m, t], m == 4];
nv3 = NeumannValue[-Fa[a, m]/\[CapitalDelta]a P[a, m, t], a == 1];
nv4 = NeumannValue[-Fa[a, m]/\[CapitalDelta]a P[a, m, t], a == 0];

(* CONSTRUCT GAUSSIAN IC*)
x0 = 1/2; (*a*)
y0 = 3/2; 
sy = 1/15;
sx =  1/100; (*a*)

h[a_, m_] = Exp[-((a - x0)^2/(2*sx) + (m - y0)^2/(2*sy) )]; 
Ni = Integrate[h[a, m], {a, 0, 1}, {m, 0, 4}];
g[a_, m_] = Exp[-((a - x0)^2/(2*sx) + (m - y0)^2/(2*sy) )]/Ni;

ufun = NDSolve[{eqn == -nv1 - nv2 - nv3 - nv4 , 
   P[a, m, 0] ==  g[a, m]}, P, {a, 0, 1}, {m, 0, 4}, {t, 0, T}, 
  AccuracyGoal -> \[Infinity], 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 1/5000}}}, 
  MaxStepFraction -> 1/200, MaxSteps -> Infinity]

imgs = Table[
   DensityPlot[ufun[a, m, t], {a, 0, 1}, {m, 0, 4}, PlotRange -> All, 
    PlotLegends -> Automatic], {t, 0, T/2, T/20}];
ListAnimate@imgs

If anyone has worked with this paper and is aware of the IC's then please let me know. I am thankful for any help with this problem.