Introduction: Gap junction nonlinear dynamics has been hypothesized to contribute to arrhythmogenicity in structurally remodeled tissue. Its implementation in large-scale 3D models remains challenging. We aim to evaluate the impact of nonlinear diffusion during normal and abnormal rhythms in a left atrial model. Methods: We generalized to 3D the Hurtado et al. mathematical framework for homogenization of conduction properties in the presence of nonlinear gap junctions. The system becomes a monodomain model in which conductivities are time-varying and depend on transjunctional potentials. Gap junction conductance was derived from the steady-state current-voltage relationship of the Vogel-Weingart model. That conductance was halved when junctional potential was > 67 mV. Neglecting the gating mechanism provided an upper bound of the effect. A bilayer interconnected cable model of the left atrium with 100 µm resolution was used. The diffusion matrix was recomputed at each time step. Normal and reentrant propagation was simulated in tissue substrates with and without random diffuse fibrosis and with reduced coupling. Simulations starting from the same initial condition were repeated with linear and nonlinear gap junctions. Results: Total activation time was imperceptibly longer with nonlinear gap junctions (145.40 vs 145.36 ms). With reduced coupling and fibrosis resulting in very slow propagation (10-15 cm/s), the cycle length of an anatomical reentry was slightly prolonged 506 ± 5 ms (nonlinear) vs 497 ± 5 ms (linear). During functional reentry in the same substrate, the cycle length was 206 ± 3 ms (nonlinear) vs 205 ± 3 ms (linear). These small differences however accumulated and progressively shifted the spiral waves. Conclusion: The effect of gap junction nonlinearity is negligible in normal conditions but can reinforce the discrete nature of propagation in slow conduction regions.