A search for new heavy resonances decaying to a pair of Higgs bosons (HH) in proton-proton collisions at a center-of-mass energy of 13 TeV is presented. Data were collected with the CMS detector at the LHC in 2016–2018, corresponding to an integrated luminosity of 138 fb−1. Resonances with a mass between 0.8 and 4.5 TeV are considered using events in which one Higgs boson decays into a bottom quark pair and the other into final states with either one or two charged leptons. Specifically, the single-lepton decay channel HH→bb¯¯¯WW∗→bb¯¯¯ℓvqq¯¯¯′ and the dilepton decay channels HH→bb¯¯¯WW∗→bb¯¯¯ℓvℓv and HH→bb¯¯¯ττ→bb¯¯¯ℓvvℓvv are examined, where ℓ in the final state corresponds to an electron or muon. The signal is extracted using a two-dimensional maximum likelihood fit of the H→bb¯¯¯ jet mass and HH invariant mass distributions. No significant excess above the standard model expectation is observed in data. Model-independent exclusion limits are placed on the product of the cross section and branching fraction for narrow spin-0 and spin-2 massive bosons decaying to HH. The results are also interpreted in the context of radion and bulk graviton production in models with a warped extra spatial dimension. The results provide the most stringent limits to date for X → HH signatures with final-state leptons and at some masses provide the most sensitive limits of all X → HH searches.