3D Diffusion with Zero Neumann Boundary Conditions

I am trying to run a time dependent diffusion model on a 3D mesh (and E. Coli bacteria) with Neumann boundary conditions equal to zero and no forcing term, f=0. I have drawn from the d2_d2D.py example problem and tried to adapt it to my problem. I have the mesh imported just fine and I think I have made the necessary adjustments for Neumann boundary conditions but when I run the model I am getting no evolution in time (i.e. my solution at t=T is the same as t=0). So I think one problem is that I can not come up with a simple example problem (that will satisfy the boundary conditions) to check if my code has an error or my problem set up. Any advice on how to go about this would be greatly appreciated. I have also put my code below if there are any glaring errors I would appreciate the feedback on that too. I am new to this forum so I am not sure how it works to share files but I can provide the .xml file with the mesh to anyone who wants to actually run the model. Thanks so much in advance for the help, I really appreciate it!!

------------------------------------
"""This demo program solves the diffusion equation

     div grad u1 = du1/dt

with only BULK diffusion of one species on the ecoli mesh with homogeneous Neumann boundary conditions and no forcing term (i.e. no reactions).
"""

from dolfin import *
import scipy.io
import numpy
#parameters["linear_algebra_backend"] = "uBLAS"

#Create mesh from Coli problem
mesh = Mesh('coli.xml')
V = FunctionSpace(mesh, "Lagrange", 1)
print mesh

#Diffusion constants for membrane and cytoplasm
Dmem = 1E-14
Dcyt = 2.5E-12

# Define initial conditions
alpha = 2; beta = 1; gamma=8;
u0 = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*x[2]*x[2] + gamma*t',
                alpha=alpha, beta=beta, gamma=gamma, t=0)
u0i = interpolate(u0, V)
plot(u0i)

T = 20 # total simulation time
dt = 0.3 # time step

# Define variational problem

# Laplace term
u = TrialFunction(V)
v = TestFunction(V)
a_K = Dcyt*inner(nabla_grad(u), nabla_grad(v))*dx

# "Mass matrix" term
a_M = u*v*dx

M = assemble(a_M)
K = assemble(a_K)
A = M + dt*K

# f term
#f = Expression('0')

# Compute solution
u = Function(V)
t = dt
while t <= T:
    #f_k = interpolate(f, V)
    #F_k = f_k.vector()
    b = M*u0i.vector() #+ dt*M*F_k
    u0.t = t
    solve(A, u.vector(), b)

    t += dt
    u0i.assign(u)

#Try to get solution evaluated at vertex points instead of cell centers
VertexSolution = numpy.zeros(mesh.num_vertices(), dtype='d')
u.compute_vertex_values(VertexSolution, mesh)

Mbulk = M.array()
Kbulk = K.array()

# export to MATLAB as a dense matrix
scipy.io.savemat("rdBULKdense.mat", {"Mbulkdense": Mbulk, "Kbulkdense": Kbulk, "VertexSolution" : VertexSolution})

file = File("BulkDiffusion.pvd")
file << u

plot(u)
interactive()
-----------------------------------------------