## 2 Computational Methods  
The simulation system consisted of a fully hydrated lipid bilayer, one PCNT or one FCNT, and 2 mol% carbon dioxide (CO\({}_{2}\)) or urea molecules. Here, mole percent is defined as the number of CO\({}_{2}\)/urea molecules divided by the sum of CO\({}_{2}\)/urea and water molecules. CO\({}_{2}\) and urea are two common metabolites, on behalf of small hydrophobic and hydrophilic molecules, and suitable for transport through CNT hollow channels. The hydrated bilayer developed by [PERSON] and [PERSON] was composed of 128 dipalmitoylphosphatidylcholine (DPPC) lipids and about 5000 water molecules. We constructed the pristine and functionalized armchair type (6,6) and (7,7) CNTs, which were about 3.6 nm in length, close to the thickness of the membrane. The FCNT was modeled by attaching 12 hydroxyl groups and 12 hydrogen atoms to the terminal carbon atoms on the pristine one.  
The force field parameters for DPPC lipids and CNTs were taken from [PERSON] et al. and [PERSON] et al., respectively [29; 30]. The topologies of CO\({}_{2}\) and urea were created by the PRODRG server, based on the GROMOS53a6 force field. However, the partial charges assigned by PRODRG were found to result in an incorrect water solubility and unrealistic partitioning between water and membrane [31]. We, therefore, resorted to Mulliken partial charges obtained after density functional theory (DFT) calculations at the B3 LYP/6-31G (d,p) level. The vdW parameters for CO\({}_{2}\) and urea were included with GROMACS itself. Water was represented by the SPC model [32]. The carbon atoms in PCNT were treated as uncharged Lennard-Jones (LJ) spheres with a cross-section of \(\sigma_{\rm cc}\) = 0.34 nm and a depth of the potential well of \(\varepsilon_{\rm cc}\) = 0.36 kJ/mol. The carbon, oxygen, and hydrogen atoms in C--O--H were assigned partial charges of +0.5\(e\), -0.8\(e\), and +0.3\(e\), the same as [23]. The bonded and non-bonded parameters of C--O--H were adopted from the GROMOS 53a6 force field [33].  
All simulations were performed under the isothermal-isobaric (NPT) ensemble by using the Gromacs package 4.5.6 [34; 35]. Periodic boundary conditions were employed in all directions. The vdW interactions were treated with a smooth cutoff at a distance of 1 nm, whereas the particle-mesh Ewald method was used to calculate the long-range electrostatic interactions [36; 37]. The temperature was kept stable at 323 K using the V-rescale thermostat, and the pressure was controlled semi-isotropically by a Berendsen barostat [38; 39]. Bond lengths within CNT/DPPC and water molecules were constrained by the LINCS and the SETTLE algorithms, which allowed a time step of 2 fs [40; 41].  
The free energy of CO\({}_{2}\) and urea molecules across the bilayer and CNT was obtained from the potential of mean force (PMF) using umbrella sampling [42]. First, we performed steered MD simulations to pull the molecule from the aqueous phase to the center of mass of the membrane or CNT [43]. Then, 35 windows were generated along the reaction coordinate (\(z\)-direction). The \(z\) coordinates of the center of mass (COM) distance between the molecule and membrane or CNT in each window differed by about 0.1 nm to ensure sufficient sampling. Each window was run for 10 ns, and the last 5 ns data were collected for sampling. In such cases, a simulation was biased by a harmonic potential with a spring of 1000 kJ mol\({}^{-1}\) nm\({}^{-2}\). Finally, the PMF profile was depicted by the weighted histogram analysis method (WHAM) [44].