This thesis delves into the potential of magnetic fusion energy, and in particular focuses on the stellarator concept. Stellarators use external coils to produce 3-dimensional (3D) magnetic fields that confine a thermonuclear plasma in a topologically toroidal volume, and thus do not, in general, require an externally driven plasma current. This is one of the main advantages of the stellarator, as the absence of strong currents in the plasma makes it intrinsically more stable than his cousin, the tokamak. It also comes at a price, since the stellarator needs to break axisymmetry to provide confinement, which is a major engineering challenge. The fusion performance in stellarators increases with β, i.e. the plasma pressure normalized by the magnetic pressure. At finite pressure however, the plasma generates additional currents, and therefore its own magnetic field, that adds up to the vacuum magnetic field generated by the external coils. The computation of the vacuum field produced by the external coils is thus not sufficient to design, optimize, operate, and interpret experimental results. Instead, it is crucial to compute the magnetohydrodynamical (MHD) equilibrium, which takes into account the non-linear contributions from the plasma. Of particular interest is the magnetic field line topology of 3D magnetic equilibria. These equilibria are in general composed of nested magnetic surfaces, magnetic islands and chaotic field lines, where the latter two topologies are, in general, detrimental to core confinement. Configurations with large regions filled with nested magnetic surfaces are thus usually sought. While it is possible to design stellarators with nested magnetic surfaces in vacuum, the plasma contribution to the total magnetic field can destroy the carefully designed magnetic surfaces at finite plasma β, thereby setting the maximum achievable β in stellarators, and ultimately limiting their performance. This thesis investigates the effect of pressure and currents generated by the plasma on the topology of magnetic field lines in MHD equilibria. Tools to compute free-boundary 3D MHD equilibria with magnetic islands and chaotic field lines are presented and extended. In particular, the Stepped Pressure Equilibrium Code (SPEC) is expanded to allow the prescription of the net toroidal current profile, and is used to perform large parameter scans to identify the equilibrium β-limits in different stellarator geometries, taking into account the effect of the bootstrap current. New measures are developed to evaluate the amount of chaotic field lines in an equilibrium, and to calculate their impact on particle transport. An analytical model is then proposed to explain the numerical results and expose the underlying scaling laws. Finally, this thesis explores the use of SIMSOPT, a python optimization framework, to optimize a configuration equilibrium β-limit. Broadly, this thesis contributes to the ongoing research on magnetic fusion reactors and the potential of nuclear fusion as a clean, safe, and abundant energy source. Specifically, it provides a better understanding of the effect of pressure on the topology of magnetic field lines in MHD equilibria, and how it impacts the performance of the stellarator. Additionally, this thesis gives insight into how optimizations can improve the performance of the stellarator and increase the equilibrium β-limit.