  The nonlinear supermatrix $\sigma $-model is widely used to understand the
physics of Anderson localization and the level statistics in noninteracting
disordered electron systems. In contrast to the general belief that the
supersymmetry method applies only to systems of noninteracting particles, we
adopt this approach to the disorder averaging in the interacting models. In
particular, we apply supersymmetry to study the Sachdev-Ye-Kitaev (SYK) model,
where the disorder averaging has so far been performed only within the replica
approach. We use a slightly modified, time-reversal invariant version of the
SYK model and perform calculations in real-time. As a demonstration of how the
supersymmetry method works, we derive saddle point equations. In the
semiclassical limit, we show that the results are in agreement with those found
using the replica technique. We also develop the formally exact superbosonized
representation of the SYK model. In the latter, the supersymmetric theory of
original fermions and their superpartner bosons is reformulated as a model of
unconstrained collective excitations. We argue that the supersymmetry
description of the model paves the way for precise calculations in SYK-like
models used in condensed matter, gravity, and high energy physics.
