Here's one way to solve the problem by explicitly modeling the interface region (note that this is explained in more detail here). We start 1st by defining some helper functions.

Mesh helper functions (edited)

These helper functions are similar to those described in the OP but have a couple of additions.

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(* Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] :=
  MeshRegion[Transpose[{data}], 
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] := 
 Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] := 
 Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] := 
 Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] := 
 Quiet@Abs@
   FindRoot[firstElmHeight[x0, xf, n, r] - fElm, {r, 1.0001, 1/fElm}, 
     Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] := 
 N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := Module[{seg},
  seg = leftSegmentGrowth[len, n, fElm];
  flipSegment[seg]
  ]
reflectRight[pts_] := With[{rt = ReflectionTransform[{1}, {Last@pts}]},
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] := 
 With[{rt = ReflectionTransform[{-1}, {First@pts}]},
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]

Build 3 region mesh

The following workflow will build a 3 material region mesh representing the left side, a thin interface region, and the right side.

(*Constants*)
l = 0.1;
hl = l/2;
int = l/1000;
hint = int/2;
hlrem = hl - hint;
nhint = 5;
nhl = 100;
intelm = hint/nhint;
(*Build right-hand side of the mesh*)
hsegint = Subdivide[0, hint, nhint];
rint = pointsToMesh@hsegint
hseg = leftSegmentGrowth[hlrem, 100, intelm/10];
rhseg = pointsToMesh@hseg
(*Extend and reflect mesh about origin*)
r = pointsToMesh@reflectLeft@extendMesh[hsegint, hseg]
(* Extract Coords from RegionProduct *)
crd = MeshCoordinates[r];
(* grab line element incidents RegionProduct mesh *)
inc = Delete[0] /@ MeshCells[r, 1];
mesh = ToElementMesh["Coordinates" -> crd, 
   "MeshElements" -> {LineElement[inc]}];
(* Association for Clearer Region Assignment *)
reg = <|"left" -> 1, "int" -> 2, "right" -> 3|>;
(*Define regions for region marker function*)
Ωleft = ImplicitRegion[-hl <= x <= -hint, {x}];
Ωint = ImplicitRegion[-hint <= x <= hint, {x}];
Ωright = ImplicitRegion[hint <= x <= hl, {x}];
(* RegionMember Function *)
rmfl = RegionMember[Ωleft];
rmfi = RegionMember[Ωint];
rmfr = RegionMember[Ωright];
regmarkerfn = 
  Piecewise[{{reg["left"], rmfl[#1]}, {reg["int"], rmfi[#1]}}, 
    reg["right"]] &;
(* Get mean coordinate of each hexa for region marker assignment *)
mean = Mean /@ GetElementCoordinates[mesh["Coordinates"], #] & /@ 
    ElementIncidents[mesh["MeshElements"]] // First;
regmarkers = regmarkerfn /@ mean;
mesh = ToElementMesh["Coordinates" -> mesh["Coordinates"], 
   "MeshElements" -> {LineElement[inc, regmarkers]}];

Set up and solve PDE system

We will use the MassTransportPDEComponentto create a PDE operator that depends on region dependent coefficients as shown below:

(*Mass transport parameters*)
dl = 10^-6;
dr = 10^-6;
kf = 0.1; kr = 1/8 kf;
K = kf/kr;
k0 = 10^4;
(*Initial condition function*)
imo[i_] := Piecewise[{{0.1, i <= -hint}, {0, i > hint}}]
iao[i_] := Piecewise[{{0, i < hint}, {0.1, i >= hint}}]
(*Define switch to activate reactions in membrane region*)
sigma = Evaluate[
   Piecewise[{{1, ElementMarker == reg["int"]}, {0, True}}]];
(*Defined region dependent diffusion coefficients*)
diffleft = 
  If[ElementMarker == reg["left"] || ElementMarker == reg["int"], dl, 
   0];
diffright = 
   If[ElementMarker == reg["right"] || ElementMarker == reg["int"], 
    dr, 0];;
(*PDE system*)
op = MassTransportPDEComponent[{{im[t, x], ia[t, x]}, t, {x}}, <|
    "DiffusionCoefficient" -> {{diffleft, 0}, {0, diffright}}, 
    "MassReactionRate" -> {{-sigma k0 (K ia[t, x] - im[t, x])/
         im[t, x], 0}, {0, 
       sigma k0 (K ia[t, x] - im[t, x])/ia[t, x]}}|>];
eqns = op == {0, 0};
{imFun, iaFun} = 
  NDSolveValue[{eqns, im[0, x] == imo[x], ia[0, x] == iao[x]}, {im, 
    ia}, x ∈ mesh, {t, 0, 60}];
imgs = Plot[{imFun[#, x], iaFun[#, x]}, x ∈ mesh, 
     PlotPoints -> 500, 
     PlotRange -> {0, Max@imFun["ValuesOnGrid"]}] & /@ 
   Subdivide[0, 60, 60];
ListAnimate@imgs

It is straightforward to change parameters and initial conditions and estimate their effect on the solution. Here, we will initialize the left-hand side at 0 concentration and give the right-hand side a much higher diffusion coefficient similar to a gas and watch the system evolve.

(*Mass transport parameters*)
dl = 10^-6;
dr = 10^-3;
kf = 0.1; kr = 1/8 kf;
K = kf/kr;
k0 = 10^4;
(*Defined region dependent diffusion coefficients*)
diffleft = 
  If[ElementMarker == reg["left"] || ElementMarker == reg["int"], dl, 
   0];
diffright = 
  If[ElementMarker == reg["right"] || ElementMarker == reg["int"], dr,
    0];
(*Initial condition function*)
imo[i_] := Piecewise[{{0, i <= -hint}, {0, i > hint}}]
iao[i_] := Piecewise[{{0, i < hint}, {0.1, i >= hint}}]
op = MassTransportPDEComponent[{{im[t, x], ia[t, x]}, t, {x}}, <|
    "DiffusionCoefficient" -> {{diffleft, 0}, {0, diffright}}, 
    "MassReactionRate" -> {{-sigma k0 (K ia[t, x] - im[t, x])/
         im[t, x], 0}, {0, 
       sigma k0 (K ia[t, x] - im[t, x])/ia[t, x]}}|>];
eqns = op == {0, 0};
{imFun, iaFun} = 
  NDSolveValue[{eqns, im[0, x] == imo[x], ia[0, x] == iao[x]}, {im, 
    ia}, x ∈ mesh, {t, 0, 60}];
imgs = Plot[{imFun[#, x], iaFun[#, x]}, x ∈ mesh, 
     PlotPoints -> 500, 
     PlotRange -> {0, Max@imFun["ValuesOnGrid"]}] & /@ 
   Subdivide[0, 60, 60];
ListAnimate@imgs

Modeling very thin membranes

@Afmo had a question in the comments about the numerical stability of modeling a 1 nm membrane on a 1 cm domain (a 7 order of magnitude difference). Now, if you're not interested in what happens within the membrane material and just want to use it as a means to model the jump condition, the membrane does not need to be 1 nm. One can simply merge the 2 fluid domains as a postprocessing step and ignore the membrane itself. However, in this case, the anisotropic meshing procedure still appears to be numerically stable over a 7 order magnitude ratio of domain size to membrane size as shown in the following workflow.

Very thin membrane mesh

To maintain numerical stability, it is generally a good practice to at least match the element thickness at the interface. An example is shown with the following workflow (please note that I adjusted the mesh helper functions to accommodate these huge aspect ratio differences):

(*Constants*)
l = 0.1;
hl = l/2;
(*7 order magnitude difference between membrane thickness and domain 
size*)
int = l/10000000;
hint = int/2;
hlrem = hl - hint;
nhint = 1;
nhl = 100;
intelm = hint/nhint;
(*Build right-hand side of the mesh*)
hsegint = Subdivide[0, hint, nhint];
rint = pointsToMesh@hsegint
hseg = leftSegmentGrowth[hlrem, nhl, intelm];
rhseg = pointsToMesh@hseg
(*Extend and reflect mesh about origin*)
r = pointsToMesh@reflectLeft@extendMesh[hsegint, hseg]
(* Extract Coords from RegionProduct *)
crd = MeshCoordinates[r];
(* grab line element incidents RegionProduct mesh *)
inc = Delete[0] /@ MeshCells[r, 1];
mesh = ToElementMesh["Coordinates" -> crd, 
   "MeshElements" -> {LineElement[inc]}];
(* Association for Clearer Region Assignment *)
reg = <|"left" -> 1, "int" -> 2, "right" -> 3|>;
(*Define regions for region marker function*)
Ωleft = ImplicitRegion[-hl <= x <= -hint, {x}];
Ωint = ImplicitRegion[-hint <= x <= hint, {x}];
Ωright = ImplicitRegion[hint <= x <= hl, {x}];
(* RegionMember Function *)
rmfl = RegionMember[Ωleft];
rmfi = RegionMember[Ωint];
rmfr = RegionMember[Ωright];
regmarkerfn = 
  Piecewise[{{reg["left"], rmfl[#1]}, {reg["int"], rmfi[#1]}}, 
    reg["right"]] &;
(* Get mean coordinate of each hexa for region marker assignment *)
mean = Mean /@ GetElementCoordinates[mesh["Coordinates"], #] & /@ 
    ElementIncidents[mesh["MeshElements"]] // First;
regmarkers = regmarkerfn /@ mean;
mesh = ToElementMesh["Coordinates" -> mesh["Coordinates"], 
   "MeshElements" -> {LineElement[inc, regmarkers]}];

Set up and solve PDE system

Set up and solve the PDE system as before:

(*Mass transport parameters*)
dl = 10^-6;
dr = 10^-6;
kf = 0.1; kr = 1/8 kf;
K = kf/kr;
k0 = 10^4;
(*Initial condition function*)
imo[i_] := Piecewise[{{0.1, i <= -hint}, {0, i > hint}}]
iao[i_] := Piecewise[{{0, i < hint}, {0.1, i >= hint}}]
(*Define switch to activate reactions in membrane region*)
sigma = Evaluate[
   Piecewise[{{1, ElementMarker == reg["int"]}, {0, True}}]];
(*Defined region dependent diffusion coefficients*)
diffleft = 
  If[ElementMarker == reg["left"] || ElementMarker == reg["int"], dl, 
   0];
diffright = 
   If[ElementMarker == reg["right"] || ElementMarker == reg["int"], 
    dr, 0];;
(*PDE system*)
op = MassTransportPDEComponent[{{im[t, x], ia[t, x]}, t, {x}}, <|
    "DiffusionCoefficient" -> {{diffleft, 0}, {0, diffright}}, 
    "MassReactionRate" -> {{-sigma k0 (K ia[t, x] - im[t, x])/
         im[t, x], 0}, {0, 
       sigma k0 (K ia[t, x] - im[t, x])/ia[t, x]}}|>];
eqns = op == {0, 0};
{imFun, iaFun} = 
  NDSolveValue[{eqns, im[0, x] == imo[x], ia[0, x] == iao[x]}, {im, 
    ia}, x ∈ mesh, {t, 0, 60}, 
   StartingStepSize -> 0.000000001];
imgs = Plot[{imFun[#, x], iaFun[#, x]}, x ∈ mesh, 
     PlotPoints -> 500, 
     PlotRange -> {0, Max@imFun["ValuesOnGrid"]}] & /@ 
   Subdivide[0, 60, 60];
ListAnimate@imgs

Obviously, one will ultimately reach a limit, but we been able to extend the anisotropic meshing to an impressive 7 order magnitude difference without evidence of great numerical instability.