The shear stress distribution in circular channels was modeled in this study using gene expression programming (GEP). 173 sets of reliable data were collected under four flow conditions for use in the training and testing stages. The effect of input variables on GEP modeling was studied and 15 different GEP models with individual, binary, tary, and quatary input combinations were investigated. The sensitivity analysis results demonstrate that dimensionless parameter y/P, where y is the transverse coordinate, and P is the wetted perimeter, is the most influential parameter with regard to the shear stress distribution in circular channels. GEP model 10, with the parameter y/P and Reynolds number (Re) as inputs, outperformed the other GEP models, with a coefficient of determination of 0.7814 for the testing data set. An equation was derived from the best GEP model and its results were compared with an artificial neural network (ANN) model and an equation based on the Shannon entropy proposed by other researchers. The GEP model, with an average RMSE of 0.0301, exhibits superior performance over the Shannon entropy-based equation, with an average RMSE of 0.1049, and the ANN model, with an average RMSE of 0.2815 for all flow depths.