  Simple assumptions made regarding electron thermodynamics often limit the
extent to which general relativistic magnetohydrodynamic (GRMHD) simulations
can be applied to observations of low-luminosity accreting black holes. We
present, implement, and test a model that self-consistently evolves an electron
entropy equation and takes into account the effects of spatially varying
electron heating and relativistic anisotropic thermal conduction along magnetic
field lines. We neglect the back-reaction of electron pressure on the dynamics
of the accretion flow. Our model is appropriate for systems accreting at $\ll
10^{-5}$ of the Eddington rate, so radiative cooling by electrons can be
neglected. It can be extended to higher accretion rates in the future by
including electron cooling and proton-electron Coulomb collisions. We present a
suite of tests showing that our method recovers the correct solution for
electron heating under a range of circumstances, including strong shocks and
driven turbulence. Our initial applications to axisymmetric simulations of
accreting black holes show that (1)~physically-motivated electron heating rates
yield electron temperature distributions significantly different from the
constant electron to proton temperature ratios assumed in previous work, with
higher electron temperatures concentrated in the coronal region between the
disc and the jet; (2)~electron thermal conduction significantly modifies the
electron temperature in the inner regions of black hole accretion flows if the
effective electron mean free path is larger than the local scale-height of the
disc (at least for the initial conditions and magnetic field configurations we
study). The methods developed in this work are important for producing more
realistic predictions for the emission from accreting black holes such as
Sagittarius A* and M87; these applications will be explored in future work.
