I want to solve two coupled partial differential equations on two dimension. There are two variables v and m. The geometry is a disk. The variable v diffuses inside the disk until it reaches the boundary and then it converts to variable m. Variable m then diffuses on the boundary, on the edge of the disk. Variable m does not exist inside the disk, it only exists on the boundary. In the diagram below you see the summary of the problem:

I use the set of equations below to define the problem:

The first equation describes the diffusion of variable v inside the disk.

The second equation describes the conversion of variable v to variable m (the term alpha*v(x,y,t)) and the diffusion of variable m on the boundary of the disk, here it is a circle.

The last equation is the boundary condition at the boundary of the disk which accounts for the conversion of variable v to variable m. On the left ∇ is the gradient operator which indicates the flux of variable v on the boundary. It will appear as the Neumann boundary condition:

NeumannValue[-1*alpha*v[x, y, t], x^2 + y^2 == 1]

Problem:

My problem is that how I should tell Mathematica that in the system of equations below (also shown above before) the first equation applies to the disk and the second equation applies to the boundary of the disk? The way I solved it below, the value of variable m is calculated on the entire of the disk which is not desired. m has value only on the boundary while it diffuses there.

Here is the code in Mathematica:

alpha = 1.0;
geometry = Disk[];

sol = NDSolveValue[{D[v[x, y, t], t] == 
     D[v[x, y, t], x, x] + D[v[x, y, t], y, y] + 
      NeumannValue[-1*alpha*v[x, y, t], x^2 + y^2 == 1],
    D[m[x, y, t], t] == 
     D[m[x, y, t], x, x] + D[m[x, y, t], y, y] + alpha*v[x, y, t], 
    m[x, y, 0] == 0, v[x, y, 0] == Exp[-((x^2 + y^2)/0.01)]}, {v, 
    m}, {x, y} \[Element] geometry, {t, 0, 10}];

v = sol[[1]];
m = sol[[2]];

ContourPlot[v[x, y, 1], {x, y} \[Element] geometry, PlotRange -> All, 
 PlotLegends -> Automatic]

ContourPlot[m[x, y, 10], {x, y} \[Element] geometry, PlotRange -> All,
  PlotLegends -> Automatic]

Adding DirichletCondition[m[x, y, t] == 0, x^2 + y^2 < 1] to enforce the value of m inside the geometry (here the disk) gives this error:

NDSolveValue::bcnop: No places were found on the boundary where x^2+y^2<1 was True, so DirichletCondition[m==0,x^2+y^2<1] will effectively be ignored.

I hope at the end I can reproduce the results of the paper below in which several proteins diffuse inside a sphere and on its surface and convert to each other on the surface. The paper is open access:

https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003396

Physical interpretation

The variable v and m represent two proteins. Protein v diffuses freely inside the cytosol (inside of the cell, here represented as a disk). Protein m is a membrane-bound protein that is it attaches to the cell's membrane (here the boundary of the disk) and only can exist as a membrane-bound protein. The protein v diffuses freely inside the disk and reaches the membrane or the boundary. There it converts to protein m with a rate that is proportional to the value of protein v on the membrane. The created membrane-bound protein m then diffuses on the membrane. Protein m cannot detach from the membrane and thus it must not exist in the cytosol (inside the disk).