I want to solve the Lokta Fisher model, a system of two differential partial equations that represents diffusion of two innovations around a City (Paris). I try with this program, but NDSolveValue don't accept to solve this system of non linear équations.

I think that I must remplace DirichetConditions by NeumannValues. But the documentation is not very clear ? Can you help me ?

(*Equation reaction diffusion temporelle *)
ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
parisradius = 0.3;
parisdiffusionvalue = 150;
carto = DiscretizeGraphics[
   CountryData["France", {"Polygon", "Mercator"}] /. 
    Polygon[x : {{{_, _} ..} ..}] :> Polygon /@ x];
paris = First@
   GeoGridPosition[
    GeoPosition[
     CityData[Entity["City", {"Paris", "IleDeFrance", "France"}]][
      "Coordinates"]], "Mercator"];
bmesh = ToBoundaryMesh[carto, AccuracyGoal -> 1, 
   "IncludePoints" -> CirclePoints[paris, parisradius, 50]];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 0.05, 
   MeshQualityGoal -> 0];
mesh["Wireframe"];

op1 = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[t, x, y]\)\) + 
   12/10*u[t, x, y] + 3/1000*u[t, x, y]*2/1000*v[t, x, y] - \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(u[t, x, y]\)\) - 
   20;
op2 = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(v[t, x, y]\)\) + 
   12/10*v[t, x, y] + 1/1000*u[t, x, y]*2/1000*v[t, x, y] - \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(v[t, x, y]\)\) - 
   10;

ini = {u[0, x, y] == 100, v[0, x, y] == 90};

ub1 = DirichletCondition[u[t, x, y] == 0, Norm[{x, y} - paris] > .6];
ub2 = DirichletCondition[u[t, x, y] == parisdiffusionvalue, 
   Norm[{x, y} - paris] < .6];
vb1 = DirichletCondition[v[t, x, y] == 0, Norm[{x, y} - paris] > .6];
vb2 = DirichletCondition[v[t, x, y] == parisdiffusionvalue, 
   Norm[{x, y} - paris] < .6];
bce = {ub1, ub2, vb1, vb2};

{usol, vsol} = 
 NDSolveValue[{op1 == 1, op2 == 1, ini, bce}, {u[t, x, y], 
   v[t, x, y]}, {t, 0, 20}, {x, y} \[Element] mesh]
usol = usol[[0]]

frames = Table[
   ContourPlot[usol[t, x, y], {x, y} \[Element] usol["ElementMesh"], 
    PlotRange -> {0, 160}], {t, 0, 20, 2}];
ListAnimate[frames, SaveDefinitions -> True]
vsol = vsol[[0]]

frames = Table[
   ContourPlot[vsol[t, x, y], {x, y} \[Element] vsol["ElementMesh"], 
    PlotRange -> {0, 160}], {t, 0, 20, 2}];
ListAnimate[frames, SaveDefinitions -> True]


NDSolveValue::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve.

Set::shape: Lists {usol,vsol} and NDSolveValue[{-20+6/5 u[t,x,y]+(3 u[t,x,y] v[t,x,y])/500000-(u^(0,0,2))[t,x,y]-(u^(0,2,0))[t,x,y]+(u^(1,0,0))[t,x,y]==1,-10+6/5 v[t,x,y]+(u[t,x,y] v[t,x,y])/500000-(v^(0,0,2))[t,x,y]-(v^(0,2,0))[t,x,y]+(v^(1,0,0))[t,x,y]==1,{u[0,x,y]==100,v[0,x,y]==90},{DirichletCondition[u[t,x,y]==0,Sqrt[Power[<<2>>]+Power[<<2>>]]>0.6],DirichletCondition[u[t,x,y]==150,Sqrt[Power[<<2>>]+Power[<<2>>]]<0.6],DirichletCondition[v[t,x,y]==0,Sqrt[Power[<<2>>]+Power[<<2>>]]>0.6],DirichletCondition[v[t,x,y]==150,Sqrt[Power[<<2>>]+Power[<<2>>]]<0.6]}},<<2>>,{x,y}\[Element]ElementMesh[<<1>>]] are not the same shape.