Solving a system of PDEs with only one Dirichlet boundary on one equation using FEM

I'm trying to solve a nonliner system of equations where there are 4 PDEs coupled through a nonlinear term j. The functions that I'm solving for are $$\phi_s$$, $$\phi_e$$, $$c_s$$, $$c_e$$, and $$j$$ using the FEniCS Python library for solving PDEs with the finite element method.

Most of the equations are defined on a 1D domain that is split into 3 equal sized subdomains, with the exception of $$c_s$$, which is defined on a 2D domain where the x dimension is partitioned into the same equal sized subdomains. The domains are coupled through the $$c_{s,e}$$ term. One extra bit of detail is that if the domains are ordered as |0|1|2|, $$\phi_s$$, $$c_s$$, and $$j$$ are only defined on domains 0 and 1.

I have a system where Neumann boundary conditions are specified on each of the equations on all of the boundaries except for $$\phi_s(x=0)$$ which has a Dirichlet boundary condition. My question is: can the functions I'm solving for have a unique solution given that only $$\phi_s(x=0)$$ has a Dirichlet boundary condition? All equations have an initial state provided at $$t=0$$, which ensure that the two time-dependent equations have unique solutions.

Unfortunately, at the moment I am unable to uniquely solve for $$\phi_s$$ in domain 2 or $$\phi_e$$ without imposing an artificial Dirichlet boundary condition (which is not technically legal in the problem space). The set of equations below can be successfully run through the COMSOL multiphysics solver without any other constraints on the problem, such as Lagrange multipliers (which are also not provided in the problem definition). I am curious what COMSOL is doing under the hood that allows it to successfully sove this problem that FEniCS is unable to do; at least in a straightforward way.

I've also tried using the Krylov solver and orthogonalizing the solution to the null space, which did allow the solver to converge on the first time step, but the solutions for $$\phi_e$$, and $$\phi_s$$ in domain 2 were mis-placed in the vertical axis (but had the correct shape).

I'm unsure of the specific details required to help me with this problem, so I can update later with any info that may be missing.

Here's the system of equations for reference:

$$\nabla\cdot\left(\sigma_{\text{eff}}\nabla\phi_s\right)=a_sFj$$

$$\frac{\partial c_s}{\partial t} = \frac{D_s}{r^2}\frac{\partial}{\partial r} \left(r^2\frac{\partial c_s}{\partial r}\right)$$

$$\nabla\cdot(\kappa_{\text{eff}}\nabla\phi_e+\kappa_{D,\text{eff}} \nabla \ln c_e)+a_sFj=0$$

$$\frac{\partial(\epsilon_ec_e)}{\partial t}=\nabla\cdot (D_{e,\text{eff}}\nabla c_e)+a_s(1-t_+^0)j$$

$$j=k_0c_e^{1-\alpha}(c_{s,\text{max}}-c_{s,e})^{1-\alpha}c_{s,e}^\alpha \left\{e^{\frac{(1-\alpha)F}{RT}\eta}-e^{-\frac{\alpha F}{RT}\eta}\right\}j$$

$$c_{s,e}=c_s(x=1, y)$$

Also, $$\eta$$ is a function of $$\phi_s$$ and $$\phi_e$$

I'd appreciate any suggestions, I've been trying to figure this out for months (I honestly don't know much about PDEs)... Thanks!