We study a nonlinear PDE-ODE system arising in bioheat modeling of localized cold exposure in which a parabolic heat equation is coupled through thermoregulatory exchange to an ODE governing core temperature, with nonlinear boundary fluxes modeling radiation, convection, and evaporation. We discuss a fully implicit finite element discretization and develop a nonlinear elliptic projection operator that accommodates the boundary coupling; we prove that this operator has uniform stability and approximation properties. These results lead next to optimal-order a priori error estimates for the backward Euler-Galerkin scheme. In the talk we focus on the construction and analysis of the projection operator, illustrate its behavior numerically, and apply it to the derivation of a priori error estimates. We illustrate with numerical results including those on multiscale modeling of the exchange terms.