The resonant substructure of D0 → π+π−π+π− decays is studied using data collected by the CLEO-c detector. An amplitude analysis is performed in order to disentangle the various intermediate state contributions. To limit the model complexity a data driven regularization procedure is applied. The prominent contributions are the decay modes D0 → a1(1260)+ π−, D0 → σ f0(1370) and D0 → ρ(770)0ρ(770)0. The broad resonances a1(1260)+, π(1300)+ and a1(1640)+ are studied in detail, including quasi-modelindependent parametrizations of their lineshapes. The mass and width of the a1(1260)+ meson are determined to be ma1(1260)+ = [1225 ± 9 (stat) ± 17 (syst) ± 10 (model)] MeV/c2 and Γa1(1260)+ = [430 ± 24 (stat) ± 25 (syst) ± 18 (model)] MeV. The amplitude model of D0 → K+K−π+π− decays obtained from CLEO II.V, CLEO III, and CLEO-c data is revisited with improved lineshape parametrizations. The largest components are the decay modes D0 → ϕ(1020)ρ(770)0, D0 → K1(1270)+K− and D0 → K(1400)+K−.The fractional CP -even content of the decay D0 → π+π−π+π− is calculated from the amplitude model to be F+4π = [72.9 ± 0.9(stat) ± 1.5(syst) ± 1.0(model)] %, consistent with that obtained from a previous model-independent measurement. For D0 → K+K−π+π− decays, the CP -even fraction is measured for the first time and found to be F+KKππ = [75.3 ± 1.8 (stat) ± 3.3 (syst) ± 3.5 (model)] %.The global decay rate asymmetries between D0 and {overline{D}}^0 decays are measured to be {mathcal{A}}_{CP}^{4uppi}=left[+0.54pm 1.04 left(mathrm{stat}right)pm 0.51 left(mathrm{syst}right)right]% and {mathcal{A}}_{CP}^{KKpi pi}=left[+1.84pm 1.74 left(mathrm{stat}right)pm 0.30 left(mathrm{syst}right)right]% . A search for CP asymmetries in the amplitude components yields no evidence for CP violation in either decay mode.