This paper focuses on the mathematical structure of the immersed boundary (IB) method, which is designed for simulating fluid-structure interaction, particularly in biological fluid dynamics. The IB method combines Eulerian and Lagrangian variables, linked by the Dirac delta function. The spatial discretization involves a fixed Cartesian mesh for Eulerian variables and a moving curvilinear mesh for Lagrangian variables. The interaction between these variables is governed by equations that use a smoothed approximation to the Dirac delta function. The paper derives the IB formulation of the equations of motion, which are based on the principle of least action. It also discusses the numerical scheme, including spatial and temporal discretization, and presents identities that ensure the conservation of mass, momentum, and energy in the discretized IB method. The paper concludes with a discussion of current and future research directions and brief applications of the IB method.This paper focuses on the mathematical structure of the immersed boundary (IB) method, which is designed for simulating fluid-structure interaction, particularly in biological fluid dynamics. The IB method combines Eulerian and Lagrangian variables, linked by the Dirac delta function. The spatial discretization involves a fixed Cartesian mesh for Eulerian variables and a moving curvilinear mesh for Lagrangian variables. The interaction between these variables is governed by equations that use a smoothed approximation to the Dirac delta function. The paper derives the IB formulation of the equations of motion, which are based on the principle of least action. It also discusses the numerical scheme, including spatial and temporal discretization, and presents identities that ensure the conservation of mass, momentum, and energy in the discretized IB method. The paper concludes with a discussion of current and future research directions and brief applications of the IB method.