We construct a heavy fermion representation for twisted bilayer graphene (TBG) systems. Two local orbitals (per spin/valley) are analytically found, which are exactly the maximally localized zero modes of the continuum Hamiltonian near the AA-stacking center. They have similar properties to the Wannier functions found in a recent study, but also have a clear interpretation as the zeroth pseudo Landau levels (ZLL) of Dirac fermions under the uniform strain field created by twisting. The electronic states of TBG can be viewed as the hybridization between these ZLL orbitals and other itinerant states, which can be obtained following the standard procedure of orthogonalized plane wave method. The ``heavy fermion'' model for TBG separates the strongly correlated components from the itinerant components and provides a solid base for the comprehensive understanding of the exotic physics in TBG.