Limit flux through subdomain boundaries

I am trying to simulate a set of partial differential equations in 2D in a material that contains multiple subdomains. I generate the geometry by reading in an image file that contains spatial information on where the edges for each subdomain are (generated from an edge recognition algorithm) and then generating edges for the PDE geometry. I attached an example image and the code for that scheme.

Ideally, I would like to set Neumann boundary conditions between these subdomains, since I can just define the faces for the geometry and set some flux condition between them. However, I know that the PDE toolbox does not allow this.

Is the only way to do this by creating a small region between each subdomain and then setting some diffusion coefficient (through the m coefficient) between these two subdomains to limit the flux reaching the other domain? This would complicate my geometry generation script and I don't really see an easy way to add an extra face between faces, since my image processing knowledge is relatively limited.

I would appreciate any ideas or suggestions on how to go about this, either from the PDE toolbox side or the image processing side.