How to define BCs on boundaries of submeshes of a parent mesh?

Hi there! I am experiencing problem defining BCs on submeshes extracted from a parent mesh.
My geometry consists of two concentric circles. I am solving one PDE for the smaller disk and another for the annulus. I successfully recovered the submeshes but having trouble defining the boundaries to my submeshes. Below is the gmsh script followed by the Python script.

gmsh script:

//+
SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 5, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 3, 0, 2*Pi};
//+
Curve Loop(1) = {2};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {1};
//+
Curve Loop(3) = {2};
//+
Plane Surface(2) = {2, 3};
//+
Physical Curve(1) = {2};
//+
Physical Curve(2) = {1};
//+
Physical Surface(5) = {1};
//+
Physical Surface(4) = {2};

Python script:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt


from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp



# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

import meshio
msh = meshio.read("example.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)


s5 = MeshView.create(cf, 5) # Extracting subspace 5
s4 = MeshView.create(cf, 4) # Extracting subspace 4

W1 = FunctionSpace(s5,"Lagrange",1) # Defining function space on s5
W2 = FunctionSpace(s4,"Lagrange",1) # Defining function space on s4

def bs5(x):  # Defining boundary for subspace s5, curve 1
  return mf == 1 
  
def b2s4(x): # Defining boundary for subspace s4, curve 2
  return mf == 2
  
# PROBLEM FOR DOMAIN 5 WITH BOUNDARY 1 
  
# Define boundary condition
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(W1, u_D, bs5)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
us = Function(W1)
solve(a == L, us, bc)

# Save solution to file in PVD format for Paraview
file = File("examplelumenvelocity.pvd");
file << us;

# PROBLEM FOR SURFACE 4 WITH BOUNDARIES 1 & 2

D, k, t = 0.7, 0.3, 3
f0, f1, beta = 5, 2, 0.3
gamma1 = (-1/D)*Constant(f0 + f1*(1-exp(-beta*(t**2)))) # Neumann condition on inner curve -dot(grad(c),n)=gamma1
gamma2 = Constant(0.0)    # Dirichlet condition on outer circle 

bc2 = DirichletBC(W2, gamma2, b2s4)

bcs=[bc2]

# Define variational problem
w = TestFunction(W2)
c = TrialFunction(W2)
a = inner(grad(c), grad(w))*dx + (((k/D)*c)*w)*dx
a -= - inner(gamma1, w)*ds_custom(bs5)
 

# Compute solution
cs= Function(W2)
solve(a == 0, cs, bcs)


# Save solution to file in PVD format for Paraview
file = File("example.pvd");
file << cs;

Error message:

*** Warning: Found no facets matching domain for boundary condition.
Invalid subdomain_id <function bs5 at 0x7fc435e3f1f0>.

How to fix this? I am using Dolfin version 2019.2.0.dev0.

I modified the BCs as

bc = DirichletBC(W1, u_D, bs5, method="pointwise") and similarly, the other one as
bc2 = DirichletBC(W2, gamma2, b2s4, method="pointwise")

Now, the error reads

ufl.log.UFLException: Invalid subdomain_id <function bs5 at 0x7f91ac0111f0>.

Does importing the parent mesh from gmsh and then extracting the subdomains mess up the physical tags or is there something wrong with my code? Any suggestion would be extremely helpful.