Atrial fibrillation (AF) is a growing clinical challenge, and computational modeling is increasingly applied to study its mechanisms and support personalized treatments. Here, we used a virtual cohort of ten bi-atrial geometries to assess how different fibrosis modeling strategies influence AF vulnerability. We developed an algorithm to generate fibrotic patterns aligned with myocardial fiber orientation and incorporated TGF-β1–mediated remodeling into the Koivumäki atrial model to simulate structural and electrical changes. A total of 200 simulations were run across four fibrosis conditions and five pacing sites. Fibrosis was modeled using four strategies: no fibrosis, decoupling only (D), decoupling with ionic remodeling (D+IR), and increased anisotropy with ionic remodeling (IA+IR). In the D model, fibrotic nodes were treated as non-conductive. D+IR combined structural decoupling with TGF-β1–induced ionic changes. IA+IR incorporated both increased anisotropy and ionic re-modeling; for this model, additional tests varied conductivity to evaluate its arrhythmogenic impact. Fibrosis levels matched Utah Stage 4, averaging 35.9% in LA and 15.8% in RA. Without fibrosis, AF was induced in only 22% of cases, underscoring the role of fibrotic remodeling. D+IR produced realistic arrhythmia rates (~46%) without tuning, while D alone led to excessively high inducibility (~60%). IA+IR required diffusion adjustments to sustain reentry (~26%). Left atrial pacing, especially near right pulmonary veins, proved more arrhythmogenic (~67%) than RA sites. These results underscore the importance of fibrosis representation and pacing location in AF modeling. Among the tested strategies, D+IR offered the best balance of physiological realism and robustness, making it a strong candidate for personalized simulations.