Starting from the Hamiltonian equations of motion of an arbitrary molecular many-body system, we first derive non-Markovian models in the form of various generalized Langevin equations (GLEs) using projection operators. The derived GLEs are integrodifferential equations for observables that are arbitrary functions of atomistic positions. We construct the projection operators to include nonlinear potentials and nonlinear memory functions in the GLEs. The primary motivation to introduce nonlinear GLEs is to move as much information as possible from the part of the GLE that ends up being modeled by a stochastic process to the deterministic part of the GLE. In this way, we ensure that one loses less information through the stochastic modeling of the exact GLE. Following this chapter, we present numerical methods to determine nonlinear memory functions from time series data. We demonstrate the numerical extraction method using a trajectory for the dihedral angle of a butane molecule in water generated by molecular dynamics simulations. From the trajectory, we calculate all previously derived GLEs using our method and compare them. For the dihedral angle dynamics of the butane molecule, we find that a position-dependent mass can lead to nonlinear memory effects. This effect can be eliminated by adjusting the mass dependence of the potential term. In the next part, we focus on the so-called approximate GLE, in which nonlinear memory effects are neglected. We discuss under which assumptions the approximate GLE emerges from a nonlinear GLE. By analytically computing the Kramers-Moyal coefficients of the approximate GLE, we show that the Fokker-Planck equation does not describe the dynamics of a non-Markovian system. We extract the friction kernel of the polypeptide Alanine9 from molecular dynamics simulations to quantify the importance of memory effects in protein folding. After parameterizing our GLE, we use the Markovian embedding method to simulate the GLE. Our GLE model very well reproduces the mean first passage times of both the folding and unfolding dynamics. The Kramers-Moyal coefficients and the mean square displacement, with pronounced anomalous diffusion, are also very well captured by the GLE. On the other hand, Markovian models based on Langevin equations with nonlinear friction cannot reproduce the dynamics in both directions with the same accuracy. From this, we conclude that consistent modeling of protein folding dynamics must take into account memory effects. The last part of the thesis is on the Markovian embedding of nonlinear GLEs. We introduce three different embedding systems that allow computationally efficient simulations of nonlinear GLEs. The first embedding system allows the simulation of nonlinear memory effects for a constant effective mass when the memory function has a nonvanishing component consisting of a delta function in time. The delta component can is not necessary for the second embedding system. We derive the second embedding from the nonlinear Zwanzig model by a perturbation expansion. The third embedding system also allows GLE simulations in the case that, in addition to a nonlinear memory function, the effective mass depends on the reaction coordinate. This embedding is not based on an approximation of the Zwanzig model and, like the first system, assumes a delta component in the memory function.