PDE in inactive form

I am trying to solve this complex PDE, however I am getting an error from Mathematica stating: The maximum derivative order of the nonlinear PDE coefficients for the Finite Element Method is larger than 1. It may help to rewrite the PDE in inactive form. I have looked at previous responses here however nothing seemed to be working when I attempted to fix it. Here is my code:

    (* ::Section::*)(*Cargo Diffusion Simulation with Reversible \
    Binding*)
ClearAll["Global`*"]
    
    
    (*1. Define parameters and functions*)
    Dcoeff = 9.24*10^-11; (*Reference diffusion coefficient (m^2/s)*)
    Xmax = 0.015;
    tmax = 10^4;
    A = Dcoeff*tmax/Xmax^2;
    Na = 6.022*10^23;
    vw = 2.99*10^-29;
    B1 = 25;
    D1 = 0.1;
    CSinitial = 1;
    CFinitial = 1/1000;
    DelU = -1;
    
    vfsb = 6.5*10^-26; (*volume factor for free,bound,and site species*)
    
    
    (*Bound cargo concentration function*)
    cb[xstar_, CF_] := (Na*cs[xstar]*CF*vw)/(Exp[DelU] + Na*CF*vw);
    
    (*Binding site concentration profile*)
    cs[xstar_] := CSinitial*(0.5 - 0.5*Tanh[B1*(xstar - D1)]);
    
    (*2. Chemical Potential Setup*)
    
    a = 2.5;
    KD = 155;
    
    (*Water volume fraction as a function of CF and xstar*)
    phiW[CF_, xstar_] := 1 - CF*vfsb - cb[xstar, CF]*vfsb - cs[xstar]*vfsb;
    
    (*Steric potential*)
    Uster[CF_, xstar_] := -(1 + a)^2*
       Log[1 - CF*vfsb - cb[xstar, CF]*vfsb - cs[xstar]*vfsb];
    
    (*Binding potential*)
    Ubind := DelU;
    
    (*Effective potential*)
    Ueff[CF_, xstar_] := Uster[CF, xstar] + Ubind;
    
    (*Spatial profile of chemical potential*)
    Mui[CF_, xstar_] := Ueff[CF, xstar]*(0.5 - 0.5*Tanh[B1*(xstar - D1)]);
    
    
    Diff2s[xstar_] := (0.7*Tanh[B1*(xstar - D1)] + 0.55 - 
        0.25*Tanh[B1*(xstar - D1)]);(*1:10*)
    
    
    (* ::Section::*)
    (*Solve PDE for C_F with nonlinear binding retardation*)
    
    With[{c = cfree[tstar, xstar]}, 
      sol = NDSolve[{(1 + (CSinitial/CFinitial)*D[cb[xstar, c], c])*
           D[c, tstar] == 
          A*D[Diff2s[xstar]*(c*D[Mui[c, xstar], xstar] + D[c, xstar]), 
             xstar] + NeumannValue[0, xstar == 0],(*Initial condition:
         same profile as before*)
         c == (0.5 - 0.5*Tanh[B1*(xstar - D1)]) /. 
          tstar -> 0,(*Boundary condition at xstar=1*)
         DirichletCondition[c == 0, xstar == 1]}, 
        cfree, {xstar, 0, 1}, {tstar, 0, 1}, 
        Method -> {"MethodOfLines", 
          "SpatialDiscretization" -> {"FiniteElement", 
            "MeshOptions" -> {"MaxCellMeasure" -> 0.005}}}]];