A cumulant expansion of the influence functional entering the path integral representation of the reduced density matrix is presented. Truncation of the expansion at the quadratic terms allows mapping of a many-body environment onto an effective harmonic model with the same force autocorrelation function, recovering the conventional Linear response limit. It is shown that the latter becomes exact in the thermodynamic limit if the bath consists of independent modes coupled to the system via coefficients that scale appropriately with the size of the bath. Low-order correction terms to the linear response approximation are found to be given by multitime correlation functions of the force exerted by the bath on the system. Within a classical approximation, such correction terms can be obtained from molecular dynamics calculations and incorporated in the path integral, which can be evaluated by means of iterative procedures. Numerical calculations on an intrinsically nonlinear bath of two-level systems confirm the above behaviors and indicate that the deviations from the linear response approximation are more significant at low temperatures where quantum mechanical effects play an important role to the bath dynamics.